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1.0  INTRODUCTION 

In  the  design  of  a  nuclear  reactor  the  critical  size  of  the  reactor  is  a 
very  important  quantity.  A  convenient  method  of  analytically  solving  the 
monoonergetic  lioltzmann  neutron  transport  equation  for  the  critical  size  is 
the  spherical  harmonics  method.  Because  the  even  order  approximations  are 
difficult  to  formulate,  most  of  the  usage  of  the  spherical  harmonics  method 
has  been  restricted  to  the  odd  order  approximations. 

Rumyantsev  (18)  was  one  of  the  first  to  suggest  that  the  use  of  the  even 
order  approximations  could  result  in  improvements  in  accuracy  over  the  odd 
order  approximations.  Mingle  (13,14)  has  found  that  for  the  disadvantage  fac- 
tor problem  the  even  and  odd  order  approximations  counterconverge  to  the 
correct  result,  i.e.,  both  the  even  and  odd  order  approximations  converge  to 
the  exact  answer  with  one  set  of  approximations  yielding  an  upper  bound  while 
the  other  set  provides  a  lower  bound.  Dawson  (4)  has  modified  the  P-  solu- 
tions to  the  Boltzmann  neutron  transport  equation  to  obtain  solutions  which 
lie  somewhere  between  the  P.  and  P  solutions.  Using  this  procedure  it  is 
found  that  the  resulting  total  neutron  flux  as  a  function  of  the  spatial  vari- 
able was  a  better  approximation  to  the  exact  distribution  than  either  the  P, 
or  P-  approximation  yielded.  In  addition  Marchuk  et,  al.  (9)  have  used  a  P- 
approximation  to  work  out  a  series  of  simple  problems.   In  most  of  these  pro- 


tho  first  order  error  in  the  critical  size  in  an  L   order  approximation  for 
a  series  of  homogeneous  slabs,  which  mathematically  predicts  that  the  critical 
thickness  determined  from  the  even  and  odd  order  approximations  will  counter- 
converge  to  the  exact  result.  There  is  thus  considerable  evidence  that  the 


even  order  approximations  will  be  of  value  in  predicting  the  critical  size  of 
a  reactor. 

In  order  to  solve  tlie  lioltzmann  neutron  transport  equation  for  a  critical 
size,  boundary  conditions  must  be  applied  at  the  outer  boundary  of  the  core 
region.  When  a  reactor  is  assumed  to  be  surrounded  by  a  vacuum  which  is  infi- 
nite in  extent,  tliese  conditions  are  termed  the  vacuum-interface  boundary 
conditions.  The  usual  vacuum-interface  boundary  conlitions  which  are  applied 
are  Mark's,  Marshak's,  and  replacement  of  the  vacuum  with  a  black  medium 
which  is  infinite  in  extent.   Pomraninc  (15,17)  has  suRRested  a  new  set  of 
boundary  conditions  developed  from  a  variational  approach  to  the  Boltzmann 
neutron  transport  equation.  These  now  boundary  conditions  have  been  found  to 
predict  the  extrapolation  distance  in  tiie  standard  Milne  problem  more  accu- 
rately than  any  of  the  boundary  conditions  currently  used  (15).  There 
appears  to  be  considerable  promise  in  the  use  of  tliese  boundary  conditions  to 
find  the  critical  size  of  a  reactor. 

The  theory,  includinR  all  four  sots  of  vacuum- interface  boundary  condi- 
tions and  computer  programs,  and  the  nature  of  the  convergence  of  the  even 
and  odd  order  spherical  harmonics  approximations  to  the  exact  critical  size 
for  various  boundary  conditions  applied  to  a  spheric.'l  reactor  is  the  subject 
matter  of  this  presentation.   In  particular  consider.able  care  is  taken  in 
developing  a  consistent  formulation  for  the  theory  of  the  even  order 
approximations. 


2.0  THEORY  OF  NEUTRON  TRANSPORT 

2.1  The  Spherical  Harmonics  Approximations 

An  extensive  development  of  the  Boltzmann  neutron  transport  equation 
can  be  found  in  many  references  (3,12,21).   In  this  work  it  will  be  suffi- 
cient to  state  that  the  general  Boltzmann  integro-differential  equation  for 
neutron  transport  can  be  written  as 
J  3/(r,E,n,t) 


-£'W(£.E,n,t)  -  {Z^(r,E,£)  +  Z^  Cr.E,£)  )^(r,E,£,t) 
S(r,E,£,t)  +  juE'jd£'^(r,E',£',t)i:^(r,E'-vE,£'-»-£)     (1) 


where:    ^(r,E,£,t)drdEd£  is  the  number  of  neutrons  in  the  differential 
volume  element  d£  about  r  having  energies  between  E  and  E  +  dE 
whose  directions  lie  in  the  solid  angle  of  phase  space  dn  about 
£,  multiplied  by  the  neutron  speed  v  (  v  =  /2E/m  ); 

!^_  (r, £'■»■£,£'•♦£) dEd£  is  the  macroscopic  cross  section  for  changing 
the  neutron  energy  from  E'  to  dE  about  E  and  the  direction  of 
motion  from  £'  to  d£  about  n,  duo  to  neutron  scattering  inter- 
actions with  stationary  nuclei; 

!:j(r,E,£)  =  dE'd£'i;^(r,E+E', £»•£')  is  the  macroscopic  scattering 
cross  section  for  any  type  of  scattering  interaction  which  causes 
a  neutron  to  be  displaced  from  its  energy,  E,  and  direction  U   at 
r;  and 

E  f_r,E,£)  is  the  macroscopic  absorption  cross  section  for  all  typos 
of  neutron  interactions  with  stationary  nuclei  which  do  not  result 
in  a  scattered  neutron. 


1    The  coordinate  system,  spherical  harmonics  and  much  of  the  nomenclature 
used  in  tliis  work  are  the  same  as  those  of  IVeinberg  and  Wigner  (21). 


The  terra  on  the  left  side  of  Eq.  (1)  is  the  time  rate  of  change  of  the 
neutron  density.  The  loss  of  neutrons  from  the  system  is  represented  by  the 
first  and  second  terms  on  the  right  side  of  Hq.  (1).  The  first  term  repre- 
sents tile  leakaj;e  from  tlie  incremental  volume  element  dr  due  to  uniform 
straightforward  motion  of  the  neutrons  whereas  the  second  term  accounts  for 
losses  of  neutrons  from  dE  about  E  and  dQ  about  n  in  energy  and  phase  space 
due  to  all  types  of  neutron-nuclei  interactions  within  the  incremental  volume 
element  dr.  The  i^ain  of  neutrons  !in  the  system  is  renresonted  by  the  third 
and  fourth  terms  on  the  ripjit  side  of  Eq.  (1).  The  third  term  accounts  for 
all  internal  and  external  neutron  sources  and  the  fourth  term  accounts  for 
the  gain  of  neutroas  resultint;  from  the  scattering  of  neutrons  into  dE  about 
E  and  d£  about  Si  from  any  other  point  in  energy  and  phase  space,  within  the 
incremental  volume;  element  dr. 

At  this  point  the  assumption  of  tlie  existence  of  only  one  monoenergetic 
thermal  neutron  grouo  is  considered.   Davison  (3)  and  otliers  (5,12,21)  have 
shown  that  this  assumption  is  valid  only  for  slightly  al)Sorbing  media  in 
regions  far  removed  from  sources  and  boundaries.  However,  even  when  these 
restrictions  do  not  apply  this  assumption  has  yielded  results  which  agree 
remarkably  well  with  experimental  results  (21).  Thus  for  simplicity  and 
without  a  great  sacrifice  in  accuracy  the  energy-independent  single  thermal 
neutron  group  is  used  throughout  the  remainder  of  this  work. 

Using  the  restrictions  that  the  individual  media  under  consideration  are 
homogeneous  and  isotropic  and  that  the  system  is  monoenergetic  and  independent 
of  time,  the  Boltzmann  equation  reduces  to 

-^•V^(r,£)  -  i^^*'-^)fi.r_,2)    *   S(r,£)  +  I d£'^(r,n")i;5 (£'-►£)  =  0  .      (2) 
Next,  defining  Z)  =  E,  ♦  I  ,  assuming  isotropic  scattering  so  that 
2j.(£'-*£J  =  ^s^'*^'   '*"''  ^"  isotropic  source  which  implies  tliat  .S(r,fi}  =  ,S(r)/4ii, 


E<\,    (2)  becomes 

-n>7/(£,£)  -  r^Cr.n)  +  ^llL  *   —  |d£'^(r,£')  =■  0  . 

In  this  treatment  only  spherical  geometry  is  considered  and  the 
coordinate  system  used  is  the  same  as  that  of  Weinberg  and  Wigner  (21).   In 
a  system  of  spherical  symmetry  the  angular  flux  distribution  depends  only 
on  two  variables:  r,  the  distance  from  the  origin,  and  u  the  cosine  of  the 
angle  6  between  the  direction  of  motion  il   and  the  extension  of  the  radius  r. 
Defining  a  new  variable  s  which  is  a  distance  laid  off  along  a,   the 
directional  derivative  in  Eq.  (3)  becomes 


(3) 


n'V/(r,n) 


^^^I'ii^  .  3£  dr  ^  a£  dp 
3s      3r  3s'   ay  3s 


W 


Figure  1.  Geometrical  interpretation  of  the  variable  s. 


Then,  as  can  be  seen  from  Fig,  1, 

ds  _  1_  ^  1_ 

3r  "  cose   u  ' 

ds  ^   r  ^     r 
38"   sine  ° 


/TV 


With  u  =  cos8, 

du  =  -sinedO  =  -/l-y-dB 

so  that 

ds  _  ds  do 
dp  "Js   du 


1-W 
and  the  directional  derivative  in  Cq.  (4)  becomes 


^Hr_,a)        ^j_^2j  D/(r,£) 


«.V/(r,£)  =  p-^_-.  --h-  ___  .  (S) 


is  also  isotropic  so  that 


S(r)    vE  *(r)    vE  |dii'/(r,S1') 

where  *(£)   =     d£'^(r,£')    is  the  total  neutron   flux  at  £.      Usini;  tlie  relation- 
ships  given   in   Eqs.    (5)    and    (6),    Hq.    (3)    becoir.es 


3^(r,£)         ^j_^2^    3^(r,£)  ^E^+i: 


Now,   d£  is  actually  d(cosO)d4i  =  dud*  so  that   by  integrating  Eq.    (7)   over  ((. 
in  the  interval  0   <  $   <  27i  and  defining 


2n 


y:(r,u)  =  I   dct./(r,£)  .  (8) 

and  tlie  mean  number  of  secondaries  per  interaction,  c,  as 


2    In  this  development  tiie  fission  neutrons  are  assumed  to  be  born  at 
thermal  energies.   In  order  to  obtain  better  results  witli  the  monoenort;otic 
assumption  the  vZ,   term  sliould  be  replaced  throughout  the  devclonment  by  a 
corresponding  terfi  representing  the  number  of  neutrons  bom  in  a  single 
fission  which  eventually  reach  tliermal  energies.  This  term  should  take  into 
account  resonance  capture,  additional  fast  fissions  and  fast  leakage.   This 
change  will  alter  the  value  of  c  but  will  not  change  any  other  quantities. 


Eq.  (7)  becomes 

In  order  to  solve  this  integro-dift'erontial  equation  the  angular  flux 
is  expanded  into  an  infinite  series  of  Legendre  polynomials,  i.e. 

^(r.g)  =  l^^MPf^M    .  (10) 

t. 

Substituting  this  expansion   into  Eq.    (9),   multiplying  the  resulting  equation 
l^y  ''m(u)<    integrating  over  the   interval   of  ortliogonality  -1    <_  u  ^■'■1.   using 
the  orthogonality  relationship 

J_^d,|.^(u)P„  (,)    =  nW  ^m   ' 
and  the   recursion   relationshijjs 


l-l"  J  — TT. TTTV   *' 


HTT 


UPJU)  = 5^;;^ 

the  spherical  harmonics  component  form  of  the  lioltzmann  transport  equation 
in  spherical  geometry  becomes 

1*1    ,1*2        d  1  ,    ,  .     I       (    l-l        d  >  .    ,  , 

*   (l.-c6^^)WJr)  =  0  .        (11) 

Equation  (11)  is  an  infinite  system  of  differential  equations  in  an 
infinite  number  of  unknowns.  Altliough  the  exact  solution  of  this  system  is 
impossible,  an  approximate  solution  can  be  found  by  assuning  that 

( ♦  -1—)  f.    .  (r)  is  negligibly  small  and  ignoring  all  higher  order  terms. 

The  resulting  approximation  will  be  called  the  P.  approximation  and  will  be 
termed  an  even  or  odd  order  approximation  according  to  the  parity  of  1., 


Increasingly  higher  orders  of  approximations  should  give  nore  accurate 
descriptions  of  the  angular  and  total  neutron  fluxes. 


2.2  The  Splierical  Geometry  Solution 

In  order  to  solve  Uq.  (11)  Weinberg  and  Wignor  (21)  propose  a  solution 
of  the  form 

^^(r)  =  (2il+l)  I   Vj(Xk)Qj(i;r/A^)  .  (12) 

Tlie  Q  (Zr/A.)'s  in  llq.  (12)  are  2/it  times  the  modified  spherical  Bessel  func- 
tions of  the  third  kind  and  are  considered  in  detail  in  Appendix  A.  Tlie 
Q  (x)'s  will  be  termed  geometricaj  functions.  Now,  substituting  the  proposed 
solution,  Eq.  (12),  into  Eq.  (11)  and  using  tlie  recursion  relationships 

,j^(x)  =<l^.^(x)  -^Q,.i{x)  , 

(^.ili)Qjx)    =Q^.^(X) 
for  tlie  Q   (x)    functions,   the   followiiii;  set  of  coupled  equations   is  obtained: 


The   recursion  relationship  of  licj.    (13)    is   similar  to  those   for  the   Legendre 
polynomials   and  the  nonsingular  part   of  the   Lui.endre  polynomials  of  the   second 


Davison   (3)    shows  tliat   tlic   (i   (X,  )'s  defined  by  the   recursion   relationship  of 
Eq.    (13)    can   be  written   as   a  certain    linear  combination   of  the    Legendre  poly- 
nomials and  the  non-singulir  part  of  the   Legendre  polynomials  of  the   second 
kind,   namely: 

8, 


"J,^"k^ 


(14) 


properly  witli  ^^i.^)  =   1.  Equation  [14)  defines  a  polynomial  in  A|^  of  dogr.-c 
4  in  which  only  alternate  powers  of  X.  appear. 

For  a  F  approximation  it  was  previously  stated  that  tlie  teria 
(-jr- +  ^)  fy^^xM   should  be  made  negligibly  small.   Referring  to  l:q.  (12)  it 
is  easily  seen  that  the  condition 

\,ii\)  -  0  CIS) 

will  satisfy  this  requirement.   Equation  (15)  defines  the  roots  A.  . 

Davison  (3)  considers  the  roots  of  Eq.  (15)  in  detail.   For  odd  order 
approximations  it  is  shown  that  the  roots  always  occur  in  pairs  such  that  one 
of  the  pair  is  the  negative  of  the  other.   For  c  <  1  all  of  the  roots  are  real 
whereas  for  c  >  1  two  of  the  roots  are  imaginary  and  the  remaining  roots  are 
real.   In  either  case  there  are  L+1  different  roots  or  (L+l)/2  pairs  of  dif- 
fejent  roots  in  an  odd  order  approximation.   For  even  order  approximations 

Davison  (3)  shows  that  there  is  always  a  zero  root.   In  addition  if 

2 
c  <  {(L+1)P^(0))   (which  is  true  for  most  nuclear  reactors)  for  a  given  even 

order  approximation,  the  remaining  L  roots  of  l;q.  (IS)  are  real  or  imaginary 

in  pairs  just  as  for  the  preceding  odd  order  approximation  roots.   There  is  thu 

a  correspondence  between  an  odd  order  approximation  and  the  succeeJin.'.  even 

order  approximation. 

The  zero  root  (X  )  in  the  even  order  approximations  muit  be  considoroti 

separately.   Following  Davison  (3)  X^   is  assumed  to  be  small  iid  tend  toward 

zero.   Then  assuming  R  is  the  interface  between  two  media  the  toiistHiii  A  is 

o 

assumed  to  be  of  the  form 
^  =  A>Q,(IR/A^)  . 

Using  this  exjiression  and  taking  the  limit  as  A  approaches  zero  of  the  tuim 
in  Eq.  (12)  corresponding  to  the  zero  root 


Um  A^(2U1)GJA^)Q|^(J:R/A^)  =  A^(2Ul)Um  G JA^)qj^(!:r/X^)/Q JER/X^) 
o  o 

.  ,A*(2)ltl)(-l)*'P  (0)  ,  r  =  R 
"  'o,  r  <  R 

which  follows  since  with  X  =  0  in  Eq.  (14),  G  (0)  =  (-Ij^T-CO).   It  is  clear 
that  this  limitinR  procedure  can  be  carried  out  only  once  witliin  a  particular 
region  and  hence  the  zero  root  term  for  even  order  approximations  is  non-zero 
at  only  one  interface  of  a  particular  region.   In  this  work  the  zero  root  term 
will  be  chosen  so  tliat  it  is  nonzero  at  tlie  right  hand  interface  of  a  region 
since  it  is  more  common  in  spherical  geometry  to  think  in  terms  of  concentric 
regions  expanding  outward  from  the  origin  rather  than  of  outer  regions  con- 
tracting inward,   llius  in  the  even  order  approximations  the  terra  due  to  the 
zero  root  is  finite  at  a  right  hand  interface  and  identically  zero  elsewhere 
in  the  medium.   This  term  introduces  a  discontinuity  in  the  even  order  moments 
at  the  interface  between  two  media. 

Finally  the  solution  of  l:q.  (II)  is,  for  odd  order  approximations 

Ltl 
^^(r)  =  (2U1)  I  \^\(\)'-li(^r/X^^)    ,  (16) 

k=l 

and  for  even  order  approximations 

f^M    =  (2U1)A^(-J)  P^(0)A(r-R)  *    (2)l*l)  I   A|,Gj(X^)Q|,(i:r/X|^)      (17) 

k~l 

where  a(r-R)  is  a  unit  pulse  function  defined  by 

'^(■^-'^  =  iJ;  r  J  R  .  (181 

These  even  order  approximation  solutions  will  be  valid  at  radii  up  to  and 

Including  the  right  hand  interfacial  boundary  since  the  discontinuity  constant, 

A  ,  which  arises  at  the  interface  is  contained  in  these  solutions.  The 
o 

solutions  in  the  next  region  will  be  valid  at  this  interface  and  at  radii  up 
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to  and  including  the  next  right  hand  interfacial  boundary  where  another 
discontinuity  will  occur,   llius  although  the  solution  of  Eq.  (11)  for  even 
order  approximations  may  be  discontinuous  it  will  be  continuously  defined 
and  finite  for  all  values  of  r. 


l:s 


2.3  The  General  Theory  of  Boundary  Conditions 

In  terms  of  the  Legendre  polynomial  expansion  of  f(.T,\i)   the  total 
neutron  flux  is  given  as 

*(r)  =.  Jdn;(r,M)  =  I  d*J  dijI^n(r)P^(M)  =  4ii^g(r)  .  (19) 

From  physical  considerations  it  is  known  that  $(0)  must  be  finite.   In  order 
to  achieve  this  condition  f   (0)  must  be  finite.  Considering  Eqs.  (16)  and 
(17)  it  is  readily  seen  that  this  can  be  accomplished  by  equating  the  Aj^'s 
corresponding  to  pairs  of  roots  of  equal  magnitude.  From  the  recursion 
relationships  for  the  G  (A.)'s,  Eq.  (13),  it  is  easily  verified  that 

G^(-V  -  (-1)\(A^)  .  (20) 

In  order  to  take  advantage  of  these  conditions  in  a  medium  which  is 
spherically  symmetric  about  the  origin  a  new  set  of  functions,  C  (x),  which 
are  the  modified  spherical  Bessel  functions  of  the  first  kind,  are  defined  as 

C^(x)  =  j{<is^M   *   (-l)*'Q^(-x)}  .  (21) 

These  new  functions  are  considered  in  detail  in  Appendix  A.  Using  these 

new  functions  the  solution  to  Eq.  (11)  in  a  region  which  is  sphericaily 

symmetric  about  the  origin  is,  for  odd  order  approximations 

(L+l)/2 
^^^(r)  =  (2U1)  I     A^Gj^(X^)Cj(Zr/X^)  ,  (22) 

k=l 

and  for  even  order  approximations 

*      >.  H^ 

fAr)   =  (2«tl)A  (-1)  Pj(0)A(r-R)  ♦  (2U1)  I   A|^G^(A^)C^(!;r/A|^)  .    (23) 

k=l 

Tlie  zero  root  term  in  even  order  approximations  does  not  Introduce  any 

difficulties  since  it  is  zero  at  r  r  0,  From  this  finiteness  of  the 


r-  14 

total  neutron  flux  at  r  "■  0  boundary  condition,  [(L'»l)/2]  constants  have 

been  eliminated. 

Next  the  boundary  conditions  which  are  to  be  applied  at  an  interface 

between  two  media  will  be  considered.  Due  to  the  discontinuity  of  the  even 

order  igoments  at  an  interfacial  boundary  in  an  even  order  approximation  it 

will  be  necessary  to  consider  the  even  and  odd  order  boundary  conditions 

separately.  In  both  even  and  odd  order  approximations  the  general  boundary 

condition  which  would  seem  to  be  the  most  desirable  is  to  demand  continuity 

of  ^(r,|i)  at  a  boundary  located  at  r  •  R.  Denoting  a  second  region  with  a  bar 

this  condition  mathematically  is 

^(R.u)  =■  T(R,Vi)  . 

or 

L  I 

I  f.mPAv)  '    I  7(R)P,(n)  .  (24) 

i-0  *    «•      4-0     * 
Multiplying  Hq.  (24)  by  P  (v)  and  integrating  over  the  period  of  orthogonality, 
-1  ^  w  ^  *1»  this  condition  reduces  to 

^j(R)  =7t(R)    t  '  0.1. •••.1-  .  (25) 

Equation  (25)  is  the  boundary  condition  to  be  applied  at  an  interfacial 
boundary  in  odd  order  approximations. 

Davison  (3)  points  out  that  in  an  even  order  approximation  only  L  condi- 

* 
tions  can  be  satisfied.  This  is  because  the  constant  A  corresponding  to  the 

zero  root  cannot  be  determined  directly  since  it  exists  only  at  an  interfacial 

boundary  between  two  media  which  is  a  point  of  discontinuity.  If  the  first 

equation  of  Eq,  (25)  is  multiplied  by  (2l+l)(-l)  P  (0)  and  subtracted  from  the 

other  L  equations  the  zero  root  tern  is  eliminated  and  the  equations 

3    Square  brackets  will  be  used  in  this  work  to  denote  the  bracket  operator, 
i.e.  [x]  -  the  largest  integer  less  than  x. 


15 


^j(R)  -  (24+l)(-l)\(0)/^(R)  =  7^(R)  -  (2Ul)(-l)''Pj^(0)7^{R) 

I  -   1,2,-..,L       C20) 
constitute  the  I,  conditions  needed  to  determine  the  remaining  L  constants. 
When  i   is  odd, P  (0)  is  zero  so  Eq.  (26)  can  be  reduced  to 
^^(R)  -  (2i*l)P^(0);^CR)  =  7(R)  -  C2I1.+  1)P^(0)7^(R) 

t  -  1,2,">.L  .     (27) 
Equation  (27)  is  the  boundary  condition  to  be  applied  at  an  interfacial  bound- 
ary between  two  media  in  even  order  approximations.  Pomraning  (16)  has  rigor- 
ously verified  the  validity  of  the  use  of  Eq.  (27)  for  even  order  approxima- 
tions using  a  variational  method  of  solution  of  the  Boltzraann  neutron  transport 
equation  with  the  Legendre  polynomial  expansion  of  the  angular  flux  as  a  trial 
function.  The  use  of  either  Eq.  (25)  for  odd  order  approximations  or  Eq.  (27) 
for  even  order  approximations  results  in  2[(L+l)/2]  conditions  at  an  inter- 
.  facial  boundary  between  two  media  which  can  be  used  to  determine  as  many 
constants. 

Since  Eq.  (27)  does  not  prescribe  continuity  of  ^  (r)  it  will  be  expected 
that  in  even  order  approximations  the  zeroth  order  moment  will  be  discontinuous 
at  an  interface  between  two  media.  This  discontinuity  is  due  to  the  zero  root 
which  is  characteristic  of  even  order  approximations.  On  the  other  hand,  in 
odd  order  approximations  continuity  of  f   (r)  at  an  interfacial  boundary  is 
assured  by  Eq.  (25).   By  Eq.  (19)  *(r)  is  proportional  to  f   ij)    so  continuity 
of  ♦(r)  is  dependent  upon  continuity  of  ^  (r) .  With  i  =  1  it  is  apparent  from 
Eqs.  (25)  and  (27)  that  ^,(r)  must  be  continuous  at  the  interface  between  two 
media  in  both  even  and  odd  order  approximations.   Letting  i  =■  0  in  Hq.  (11) 
it  is  seen  that  ^,(r)  is  proportional  to  the  derivative  of  f   (r),  und  tlierefore 
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boundary  requires  that  the  neutron  current  which  is  proportional  to  the 
derivative  of  the  total  neutron  flux  bo  continuous  across  the  boundary. 
Suniarizing  the  conditions  listed  above,  t(r)  is  continuous  across  an  inter- 
facial  boundary  between  two  media  for  odd  order  approxiaations  but  discontin- 
uous across  the  boundary  tor  even  order  approximations  whereas  the  neutron 
current  is  continuous  across  a  boundary  in  both  even  and  odd  order  approxiaa- 
tions. Tlieso  conditions  can  be  seen  graphically  for  the  P.  and  I'  approxima- 
tions in  Fig.  2.  I 


Total 
flux, 
♦  {r) 


radius,  1' 


Figure  2.  The  spatial  dependence  of  the  total  neutron  flux  for  a 
two  region  problea  with  the  interfacial  boundary  at  R. 


Fron  Fig.  2  it  is  readily  apparent  that  the  even  and  odd  order  approxiin 
tions  tend  to  approach  the  exact  solution  for  the  total  neutron  flux  in  a 
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counterconvergent  manner  i.e.  the  odd  order  approximations  approach  the  exact 
solution  fron  below  in  the  central  region  while  the  even  order  anproximations 
approach  the  exact  solution  from  above.  Thus  it  is  expected  that  if  the 
central  region  was  that  of  a  reactor  core  and  the  peripheral  region  an  infi- 
nite reflector  the  computed  critical  radii  from  the  even  and  odd  order  approx- 
imations will  be  counterconvergent.  Davison  (2)  has  derived  a  formula  which 
predicts  this  type  of  counterconvergence  in  slab  geometry  for  a  series  of 
slabs.  For  values  of  c  close  to  unity  his  formula  is 

R  -  R  .  i^i-^)''*l)f    |i  ,  0{inl2Li31)}  (28) 

''    96(L+3/2)'^  '' 

where  K  is  the  exact  critical  radius  and  R.  is  the  critical  radius  found  using 
the  L   order  approximation.   Although  this  equation  was  originally  derived 
for  slab  geometry  cases  the  same  counterconvergent  trend  should  appear  in 
spherical  geometry  for  interfaces  at  which  the  curvature  of  the  interface  is 
small  with  respect  to  the  critical  radius  of  the  central  region.   A  counter- 
convergence  of  the  even  and  odd  order  approximations  will  be  expected  whenever 
an  interface  between  two  media  is  the  dominant  boundary,   Davison  (3)  has 
shown  that,  at  least  for  slab  geometry,  the  vacuum- interface  boundary  condi- 
tions attributed  to  Mark  are  equivalent  to  replacing  the  vacuum  with  an  infi- 
nite black  medium.   Since  the  replacement  of  the  vacuum  by  an  infinite  black 
medium  introduces  an  interfacial  boundary  between  two  media,  counterconvergence 
of  the  even  and  odd  order  approximations  will  also  be  expected  for  Mark's 
vacuum-interface  boundary  conditions. 

Next,  the  boundary  conditions  at  a  right  hand  vacuum  interface  will  be 
considered.   At  a  vacuum  interface  at  r  »  R  the  exact  boundary  condition  is 

/(R,m)  =0   for   M  <  0  .  (29) 

This,  however,  constitutes  an  infinite  number  of  conditions  which  cannot 


all  be  exactly  satisfied  in  an  approximation  of  finite  order.   In  an  L   order 
approximation  in  a  central  renion  only  [(L+l)/2]  conditions  can  be  satisfied 


condition,  Eq.  (29),  is  commonly  reduced  to  [(L+l)/2]  conditions  in  four  ways. 

First,  it  is  physically  plausible  tliat  tlie  total  number  of  neutrons 
entering  from  tlie  vacuum  stiould  be  zero,  'i'his  condition  can  be  written  mathe- 
matically as 

^(H,M)udu  =  0  .  (30) 

'-1 

Since  u  is  just  f'.(ii),  for  an  L   order  approximation  a  set  of  boundary  condi- 
tions which  includes  V.q.    (30)  is 

or,  what  is  equivalent 

i: 

These  vacuum-interface  boundary  conditions  are  l^nown  as  the  Marshal;  boundary 
conditions.  Choosing  the  form  in  Uq.  (32),  ex|ianding  ^(R,u)  by  Eq.  (10),  and 
interchanging  the  order  of  integration  and  summ.ition  Hq.  (32)  becomes 


j:(R,u)P,^_^(u)di,   =   0  j    =    1,2,---,[(L>-1)/2|  (31) 


^(K,M)M'J"'dM    =   0  )    =.    l,2,---,[iL*l)/;:]    .  (32) 


I  .0  , 

I  faW]      P,(m)m'^      dM  =   0  j    =    l,2,-.-,|(Ltl)/2)    .  (33) 

1=0   ^       J-1    * 

Second,   Mark    (11)    hai   [iroposed   tliat   liq.    (21))    be   satisfied  exactly   for 

[(L+l)/2]   values  of  p.      Ii  particular  he  chose  to  use  the  condition  that 

^(K,Mj)    =0  j    =    l,2,-..,|(Ul)/2]  (34) 

where  the    m's   are   defined  by 

Vl^l^j^    =0      ,      p.    <   0    .  (35) 
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From  the  Legendre  polynomial  expansion,  Hq.  (10),  it  is  readily  apparent  that 

Mark's  boundary  conditions  impose  the  condition  that  the  term  Pi  ^i  (i')^l*i('*) 

be  zero  and  thus  that  the  L*l  order  term  and  assumably  all  higher  order 

terms  in  the  Legendre  polynomial  expansion  make  no  significant  contribution 

to  the  angular  flux  at  the  vacuum  interface.  Expanding  ^(R.iJi)  by  Eq.  (10), 

Eq.  (34)  becomes 

L 

I   f AR)P Au.)    =0     j  =  1,2....,[(U1)/21  .  (36) 


J.=.0 


i.^"'-  l^-j 


Third,  the  vacuum  can  be  replaced  by  a  fictitious  infinite  black  medium. 
For  slab  geometry  Mark  (11)  has  shown  that  this  condition  is  equivalent  to 
Eq.  (36).  However,  it  is  impossible  to  prove  equivalence  of  those  two  sets 
of  boundary  conditions  for  spherical  geometry  due  to  the  curvature  of  bound- 
aries and  inherent  boundedness  in  one-dimensional  spherical  geometry.  Physi- 
cally this  difference  occurs  in  spherical  geometry  because  the  infinite  black 
medium  reflector  boundary  condition  takes  the  curvature  at  the  boundary  into 
account  in  a  much  different  manner  than  do  Mark's  vacuum-interface  boundary 
conditions.   In  general  the  infinite  black  nedium  surrounding  a  center  region 
constitutes  a  two  region  problem  in  which  c  =  0  in  the  infinite  outer  or 
surrounding  medium.   As  before  it  is  to  be  noted  that,  although  the  problem 
here  considers  only  spherical  geometry  and  Davison's  proof  of  counterconvergence 
is  for  slab  geometry,  the  same  counterconvergence  trend  should  be  found.  This 
should  be  true  for  Mark's  boundary  conditions  as  well  as  for  the  fictitious 
infinite  black  reflector  boundary  conditions. 

Fourth,  Pomraning  (15)  has  used  a  variational  formulation  of  tlie  Boltzmann 
equation  and  a  Legendre  polynomial  expansion  of  the  angular  flux  as  a  trial 
function  to  determine  a  set  of  mathematically  consistent  boundary  conditions. 
By  setting  the  first  variation  of  the  boundary  terms  to  zero  I'omraning  found 
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the  boundary  conditions  which  result  from  the  variational  calculus  to  be: 

^(R,U)«^(R.-U)pd(i  =  0  ,  (37) 

''-1 

These  boundary  conditions,  called  variational  boundary  conditions,  are 

considered  in  detail  in  Appendix  B.   In  general  they  can  be  put  into  the  form 


L 

I   /JR)b,.  =0     j  =  1.2,---,[(Ltl)/2)  .  (38) 

i=0 


r^l-^^^j 


Considering  Eqs.  (33),  (35)  and  (38)  it  is  readily  seen  that  each  of  the 
equations  can  be  forced  into  the  form  of  Eq.  (38).   For  Marshak's  vacuum- inter- 
face boundary  conditions  tliis  is  accomplished  by  setting 

b^j  =  I  P^(u)M^''^dii  (39) 

while  for  Mark's  vacuum-interface  boundary  conditions 

\i    -   P,(Uj)  .  (40) 

Since  the  variational  boundary  conditions  are  developed  directly  from 
the  mathematics  instead  of  from  adaptations  of  physical  considerations  to  the 
mathematics  as  Marshak's  and  Mark's  boundary  conditions  are,  the  variational 
boundary  conditions  should  give  the  best  results  for  a  given  order  of  approxi- 
mation.  For  low  orders  of  approximations  the  Marshak  boundary  conditions 
should  give  better  results  than  Mark's  or  the  infinite  black  reflector  bound- 
ary conditions  for  the  same  order  of  approximation  jince  the  physically 
plausible  condition  of  no  return  current  of  neutrons  from  the  vacuum  is  auto- 
matically included  in  Marshak's  boundary  conditions  whereas  it  is  not  included 
in  the  other  boundary  conditions. 
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2,4  Formulation  of  Boundary  Conditions 

In  this  section  three  liypothetical  reactor  systems  will  be  considered. 
For  each  reactor  system  the  appropriate  boundary  conditions  will  be  imposed. 
Then  the  equations  representing  the  appropriate  boundary  conditions  will  be 
forced  into  a  common  matrix  equation  form.   Finally  the  solution  of  the  common 
matrix  equation  form  will  be  examined. 

A  hypothetical  reactor  consisting  of  a  spherical  central  core  region  with 
c  >  1  surrounded  by  a  vacuum  infinite  in  extent  will  be  referred  to  in  this 
work  as  a  bare  spherical  reactor.  Choosing  the  coordinate  system  so  that  the 
center  of  the  bare  spherical  reactor  is  at  the  origin,  liqs.  (22)  and  (23)  are 
the  appropriate  solutions  of  Kq.  (11)  for  r  <  R.  The  vacuum  is  not  capable  of 
supporting  a  neutron  flux  and  therefore  cannot  be  considered  a  medium.  Since 
an  interfacial  boundary  between  two  media  does  not  exist  in  the  system  consid- 
ered here  the  constant  A  will  be  assumed  to  be  zero.  The  consequences  of 
this  assumption  will  be  examined  later  from  tin;  results  of  numerical  computa- 
tions bused  on  this  assumption.   This  assumption  allows  iiqs.  (22)  and  (23)  to 
be  simplified  to 

[(Ul)/2] 
^^(r)  =  (2..1)  J^  \^,(\i\i^^/\)    ■ 

Ajjplying  the   general  vacuim-interface  boundary  condition,   Eq.    (38),   the  set 
of   t(L+l)/2]    linear  e 
spherical   reactor  is 


l(Ul)/2]      L 

I       A     I    (2lt*l)G    (X  )C    (r.R/A    )b        =   0 
k=l        "i^O  IKK.  K      nj 


k=l       "i= 

3  =  l,2,...,[(Ul)/2]  .       (41) 
Or,  defining  the  [(L+l)/2]  by  one  column  matrix  X  with  elements 
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and  the  [(L+l)/2]  by  l(L+l)/2]  square  matrix  T  with  elements  t..  such  that 

Eq.  (41)  can  be  rewritten  in  matrix  form  as 

TX  =  0  .  (44) 

Equation  (44)  is  the  general  matrix  form  in  which  all  of  the  cases  to  be  con- 
sidered in  this  work  will  be  written.  For  a  bare  spherical  reactor  Eqs,  (42) 
and  (43)  define  the  elements  of  the  X  and  T  matrices  involved  in  Eq.  (44). 

A  hypothetical  reactor  consisting  of  a  spherical  central  core  region  with 
c  >  1  surrounded  by  a  black  medium  (c  =  o)  infinite  in  extent  will  be  referred 
to  in  this  work  as  a  bare  spherical  reactor  with  an  Infinite  black  reflector. 
Choosing  the  coordinate  system  so  that  the  center  of  the  reactor  is  at  the 
origin,  Eqs.  (22)  and  (23)  are  the  solutions  of  Eq.  (11)  for  r  <^  R.   For  r  >^  K 
it  is  known  from  physical  considerations  that  as  r  approaches  infinity  the 
total  neutron  flux  approaches  zero  i.e. 

4im  ♦(r)  =  Him  4ii/  (r)  =  0  .  (45) 

r-*«>       r-K* 

Thus  all  of  the  A.  's  corresponding  to  positive  roots  must  be  identically  zero. 

Since  for  the  infinite  black  medium  region  no  right  hand  boundary  exists  the 

T  constant  will  be  zero,  i.e.  it  will  never  arise  at  least  not  for  finite  r. 
o  ' 

Thus  for  r  >^  K  eqs.  (16)  and  (17)  reduce  to 

[(L+l)/2]    _        _   ■ 
T^M   =  (2U1)(-1)''    I       \^^iX^)q^l-TT/X^)  (46) 

where  X.  is  assumed  to  be  positive,  for  both  even  and  odd  order  approximations. 
In  order  to  solve  for  the  constants  (A.  's  and  A,  's)  in  odd  order  approximations 
Eq.  (25)  is  applied  at  the  interface  at  r  =  R  between  the  reactor  and  the  black 
medium  so  that  for  a  bare  spherical  reactor  with  an  infinite  black  reflector. 


[(L*l)/21 

(2i+i)  I     A^r,^(\|^)c^(i;R/x^) 

.[CL+l)/2] 

-  (24+l)C-l)     I         \'^4(\)QtC-^'*/*k^  °  "    £  =  0.1,-.-,L  .    (47) 

k=l 

For  even  order  approximations  the  corresponding  condition,  Bq.  (27),  is 
applied  at  the  interface  so  that 

((L+l)/2] 
(2i+l)    I        {GJXj^)C^(ER/A^)  -  Pj(0)C^(ER/X^))A^ 
k«l 

[(Ul)/21  _  ^   _ 

-  (2i+l)   I   {(-1)\(X^)Q^(-WX^)  -  PJO)Q^(-IR/T^)}A^  -  0 

I  -   1,2,.--,L   .  (48) 

Or,  defining  the  2[(L+l)/2]  by  one  column  matrix  X  with  elements 
x^  =  A^  for  k   <   [(Ul)/2]    , 

Xy.  =  I.  for   [(L*l)/2]    <  k   <  2[(L+l)/2]  (49) 

where 

j   =  k   -    [(Ul)/2] 
and  the  2[(L+l)/2]   by  2[(L+l)/2]    square  matrix  T  with  elements  t^^  such  that, 
for  odd  order  approximations 

t^^.   =    (2il*l)Gj^(X^)C|^(ER/\)  for  k    <_  [(Ul)/2] 

I  =   0,l.v,L 
t^^  =   (2Ul)(-l)'^j(T.)Q^(-Ik/T.)  for   [(L*l)/2]    <  k   <  2[(Ul)/2] 

i.  =  0,1, .••,L  (SO) 

Where 

j    =  k   -    [(Ltl)/2] 
and   for  even  order  approximations 
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for  k  _<  [(L*l)/2] 
4  =  1,2,"-,L 

t^i^  =  (2oi)(-i)'^jT.)Q^c-rR/r.)  -  (2i*i)p^(o)Q^(-rR/x.) 

for  [(Ul)/2]  <  k  _<  2[(L*l)/2) 

i   =  1,2,---,L  (SI) 


where 


j  =  k  -  [(Ul)/2] 
Eqs.  (47)  and  (48)  can  be  written  in  the  same  matrix  form  as  Bq.  (44).   For  a 
bare  spherical  reactor  with  an  infinite  black  reflector  Eqs.  (49)  and  (50)  or 
(51)  define  the  elements  of  the  X  and  T  matrices  for  Eq.  (44). 

A  hypothetical  reactor  consisting  of  a  spherical  central  core  region  witli 
c  >  1  surrounded  by  a  concentric  spherical  shell  of  thickness  R  in  which  c  <  I 
and  immersed  in  a  vacuum  infinite  in  extent  will  be  referred  to  as  a  reflected 
reactor.   Choosing  the  coordinate  system  so  that  the  center  of  the  reflected 
reactor  coincides  with  the  origin,  liqs.  (22)  and  (23)  are  the  solutions  of 
Eq.  (11)  for  r  <  R.  As  in  the  case  of  the  bare  spherical  reactor  the  zero 
root  constant  in  even  order  approximations  associated  with  the  vacuum  boundary 
will  bo  assumed  to  vanish.  Since  the  zero  root  terms  in  even  order  approxima- 
tions are  taken  to  be  nonzero  only  at  right  hand  boundaries,  the  zero  root  terra 
will  vanish  entirely  from  the  reflector  region  solutions.   Applying  the  inter- 
facial  boundary  condition,  for  odd  order  approximations,  Eq.  (25),  and  Eq.  (27) 
for  even  order  approximations,  at  r  =  R  and  the  general  vacuum- interface  bound- 
ary condition,  Eq.  (38),  at  r  =  R  +  R  the  equations  determining  the  A, 's  and 
X  's  for  a  reflected  reactor  are,  for  odd  order  approximations 


[(Ul)/2] 

2[(L*1)/21_ 

-  (21+1)    I        \i^i(\)'J;t(!:i«/\)  =  0 
k=l 

I   =  0,1,...,L       (52) 

and  for  even  order  approximntioas 

[(Ul)/2) 
(2U1)    I        {i^^(,\^iL^l,lH/X^)    -    1>J0)(:^(X1</X|^))A^ 

[U.*l)/2| 

-  (2)tti)   I      c,^(A|.)q^(!:r/X^.)  -  i';t(i')Qo(!:R/X|^)}I^  =  o 

k=l 

I  =    1,2,---,L  (53) 

together  witli 

2[(U1)/21_     I, 

I        A    I  (2t*i)(T  (X  )()  (!;(ici(  )/T  )F      =  0 

k=l  '^«.=  ()  t      k       «  OKI., 

j    =    l,2,..-,[(Ul)/2]    .  (54) 

Now  define  the   3[(1.*1)/2|    by   omj   colomii   matrix   X  witii  elements 
X|.    =   Aj^  for  k   <_   l(Ltl)/2|, 

x.    =  J.  for    [(Ul)/21    <   k    <   3|(Ulj/.;|  (55) 

^  }  ~ 

whore 

j    =   k   -    [(Ul)/2) 

and   tlie    3[(L*l)/2]    square   matrix  T  'Aith  elements   1 1^,  ■  such  that,    for  odd  order 

approximations 

l*l)G^(X^)L^iUi/>.^)  for 

Jl  =   0 , 1 ,  ■  • .  ,  L 

t^^  =    (2;ttl)TTj^(T  jQj^(rR/I.)  for    [(l.*l)/2|    <   k    <   3((Ltl)/21 

I   =   0,1, •••,L  (56) 
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wliere 

j   =  k   -    [(Ul)/2] 
and  for  even  order  approximations 

for  k  jc   [(L+l)/2] 
t  »   1,2,...,L 
t^^  =  (2i+l)(Ij^{Tj)q^(a/X.)  -  (2ul)P^(0)Q^(in</X.) 

for   [(Ul)/2]    <  k   <^3[(L*l)/2] 
I  =   1,2,...,L 
where 

j   =  k  -    [(Ltl)/2]    . 
Then  using  Eqs.    (56)   and   (57)   together  with 
t^^  =0  for  k   <   I(L+l)/2] 

2[(L+l)/2]    <  t   <^  3[(Ul)/2] 


(57) 


(58) 


where 


t,,  -     I    (2™.l)?jX.)(^(r(R.R^)/X.)F^ 


for    [Cl.+  l)/2)    <  k   <  3[(Ul)/2] 
2l(Ltl)/2]    <   I   <  3[(Ul)/2] 


(59) 


j    =  k   -    [(Ul)/2]    , 

n   =   t  -  2[(L+l)/2] 
Eqs,    (52),    (53)   and   (54)   can  be  written  in  the  same  matrix  form  as  Eq.    (44). 
For  a  reflected  reactor  Eqs.    (55),    (55)   or   (57),    (58)   and   (59)   define  the 
elements  of  the  X  and  T  matrices   in  Eq.    (44). 
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Since  all   three  of  the   reactor  systems   considered  have  been  put   into  tlie 
form  of  Eq.    (44),   solution   of  this  type  of  equation  will   be  considered  in 
detail.      I'rom  the  theory  of   linear  alRebra   (22)    the  column  matrix  X  can  liave 
a  non-trivial   solution   only  if  the  determinant  of  T  vanishes.      Hxamining  the 
various  terms  wliich  make  up  T  it   is  readily  apparent  that   It  is  the  undeter- 
mined  free  quantity  which  must   be  varied  until 

|T|    =  0   .  (00) 

Since  the  t^,(x)    functions  can  be   ijeriodic  an   infinite  number  of  values  of  R 
will  satisfy  liq.    (60).      In  order   for  tlie  total  neutron   flux  to  be  non-negative 
(a  physical  necessity)   the   first   or  smallest   positive  value  of  R  which  satis- 
fies liq.    (60)    is   tne  required  value  and   is   commonly  known  as  the  "critical" 
r£idius  of  the  spherical   reactor   in  question,      liquation   (60)    is  known  as  tlio 
criticality  equation. 

After  a  "critical"   radius   lias   been    found   the   constants   composing  the 
column  matrix  can   be   determined   to  within   a  normalization    factor.      If  the    first 
column   is  transferred  to  the  other  side  of  the  equality  and  any  one  of  the 
rows  deleted,  the   resulting  matrix  equation   (one   less   in  rank)   may  be  solved 
for  the  remaining  A^ ' s   and  ^u'^   in  terras  of  A   .     Using  this  procedure  A.    is 
defined  to  be  equal   to  one. 
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3.0   UISCUSSION  AMD   WHSULTS 


3.1  General  Discussion 


Tlie  basic  purpose  of  this  study  is  to  examine  the   rate  of  convergence  of 
the   splierical   harmonics   method  for  botli  bare  and  reflected  spheres  usinj; 
various  vacuum-interface  boundary  conditions.     Since  an  adequate  treatment   of 
the  problem  requires  a  great  number  of  numerical   computations,   a  computer  pro- 
gram is  used  to  carry  out   the  numerical  work.     The  program  is  written  in  a 
general   manner  so  that   only  one  program  is   required   for  all  of  the  cases  to  be 
considered.     The  program  is  written   for  the   IBM  1410  in  the   FORTRAN   language 
and  is  compiled  in  a   fourteen  digit   floating  point  mantissa   length.      In  the 
sections  of  the  program  which  require   iteration  a  relative  accuracy  of  lO"   ^ 
is  employed.     Since  only  six  digits  are   retained   in  the   final   results,   this 
accuracy  is  deemed  satisfactory   for   the  size  of  matrices  encountered  in  the 
program.      A  complete  description  of  tlie  computer  program  is   given  in  Appendix  C. 

Before   considering   some  practical   cases   the   equations  will   be   analyzed   for 
the   important  unspecified  constants.      It   is   readily  apparent   from  the  tlieory 
for  a  bare   sphere  tliat   c   and  £  are  the  only  unspecified  constants.     The  con- 
stant  i;  always  occurs   in  a  product   with  the   radial   variable   r  so   if  radii   are 
measured  and   reported  iji  mean   free  paths,   only  the  mean  number  of  secondaries 
per  collision  c  needs   to  be   specified   for  a  bare   spherical   reactor.      I'or 
reflected  spheres   c,    E,    c",  f,    and  R^   are   the  unspecified  constants.      Tim   con- 
stants  t  and  E  always  occur  in  products  with   r  so  if  radii  are  measured   in 
mean   free  paths,   only  Z/J,   c,   c"  and  R     are  the  unspecified  constants   for  a 
reflected  reactor. 

In  order  to  fully  appreciate  the  magnitude  of  certain  transport   tlieory 
constants  they  will   be  compared  witli  the  more   familiar  diffusion  theory 
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constants.   In  particular  tlie  thickness  of  a  reflector  in  diffusion  lengths  is 
a  very  important  property  in  diffusion  theory.   Glasstone  (6)  shows  that  if  a 
reflector  is  1.5  to  2.0   diffusion  lengths  in  tliickness  it  may  be  considered 
essentially  infinite.   From  diffusion  theory  the  diffusion  length,  <'    ,   is 

<'^  =  /uTT  =  i/ZxTT  (61) 


sion  tlieory  the   root   of 

I 


-2 


1/ 


/^)-c) 


or,    since    for  a  non-multiplyini;  mediura  c   =   Ijll, 


(62) 

With   U     being  tlie  thickness  ot   the   reflector  in  mean   free  paths 

ZV.  /J,    =  1!  Ay  I  =  7k  (63) 

0      1  0        a  o  '■      ' 

is   the  thickness  of  tlie   reflector   Iei  diffusion    lengths.     As  tlie  order  of 
approximation   increases   the   jirincipal   or   largest    root    of  bq.    (15)    clianges 
slightly  but  the  quantity  Tu  /Y    st 
the  reflector  in  diffusion    lengths. 

At   this   point    furthei-   consideration  will   be   given   to   Hq.    (28)    which  pre- 
dicts counterconvergence   of  the  oven  and  odd  order  ap|)ro>.imations   for  inter- 
facial   boundaries,    infinite  black   media   reflectors   and   Mark's   vacuum-interface 
boundary  conditions.      If  the  0( —     j        }   terms   in  liq.    (28)    arc   ignored   it   is 
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possible  to  find  a  pair  of  weighting  factors  which  would  determine  a  weighted 
average  critical  radius  that  will  approxiraate  the  exact  critical  radius.  With 
L  odd  this  combination  of  corresponding  even  and  odd  order  results  can  be 
written  mathematically  as 

where 

R-K,  , 
w  =    h*'- 


L+1  "  It, -R,  ,  • 
L   Ltl 

The  values  of  the  weighting  factors  w.  and  w   .  are  worked  out  readily  by 

using  Eq.  (28).   Ignoring  the  higher  order  terms  these  factors  are  found  to  be 

,2 


and 


2(2L*3) 


''   2(2L*3)^  ♦  (2L+5)' 


(2L+5)' 

""75 


^*^        2(2L+3)^  ♦  (2L+5)^ 


(65) 


A  list  of  the  values  of  these  weighting  factors  for  approximations  up  to  order 

nine  may  be  found  in  Table  I. 

Table  I 

Weighting  Factors  for  Eq.  (64) 

L  w,  w,  , 

L  Ltl 


1 

.505051 

.494949 

3 

.572438 

.427562 

5 

.600355 

. 399645 

7 

.615548 

.384452 

9 

.62508y 

.374911 
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For  very  large  orders  of  approximations 
Him  w  =  2/3  , 

Him  w    =  1/3  . 

Thus  for  very  large  orders  of  approximations  it  will  be  expected  that  the  odd 
order  approximations  will  be  about  twice  as  close  to  the  exact  answer  as  their 
corresponding  even  order  approximations.   Although  this  criterion  was  origi- 
nally developed  for  slab  geometry  it  should  apply  to  the  spherical  geometry 
cases  considered  in  tliis  work  when  the  critical  radius  is  large  enough  so  that 
the  curvature  of  the  boundary  contributes  a  small  effect. 
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3.2  Results  of  Numerical  Calculations 
for  Uare  Spherical  Keactors 

In  order  to  examine  the  rate  of  convergence  of  the  spherical  harmonics 
approximations  for  various  boundary  conditions  imposed  at  the  surface  of  a 
bare  sphere  three  cases  will  be  considered.   First  a  bare  spherical  reactor 
with  c  =  1.05  in  the  core  region  is  utilized.   In  this  case  the  value  of  c  is 
probably  too  large  for  diffusion  theory  to  be  exact  but  not  so  large  that  the 
assumption  of  the  existence  of  a  single  thermal  neutron  group  will  no  longer 
be  valid.   This  basic  intermediate  case  will  be  considered  in  detail  and  the 
results  from  other  cases  will  be  compared  to  the  results  of  this  basic  case. 
'Die  second  case  to  be  considered  will  be  that  of  a  bare  spherical  reactor  in 
which  c  =  1.02.   For  this  case  the  value  of  c  is  probably  close  enough  to 
unity  so  that  diffusion  theory  should  give  very  good  accuracy.   As  a  third 
case,  the  extreme  of  a  large  value  of  c,  nainely  1.40,  will  be  considered.   In 
order  to  adequately  describe  this  bare  spherical  reactor  case  a  multigroup 
transport  theory  would  bo  needed  since  the  fast  leakage  would  be  so  great  that 
the  single  thermal  energy  neutron  group  assumption  would  not  be  valid.  How- 
ever for  a  limiting  st^dpoint,  the  results  i:or  this  large  value  of  c  are 
needed  in  this  work. 

Table  II  is  a  compilation  of  the  computed  results  for  a  bare  spherical 
reactor  with  c  =  1.05.   The  exact  critical  radii  given  in  this  table  and  other 
tables  are  taken  from  results  computed  by  Carlson  and  Bell  (1)  using  the 
extrapolated  end  point  method.   The  per  cent  error  column  is  the  percentage 
error  in  an  L   order  approximation  which  is  computed  as 

R-R 
per  cent  error  =  — j, —  x  lOU  . 


The  column  labeled  "inward  flux"  represents  the  integral  of  the  computed  angular 
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Table  II.   Calculated  results  for  a  bare  spherical  reactor  in  whicli  c 
is  1.05.  The  exact  critical  radius  is  7.2772  mean  free  paths. 


Marshak  vacuum-interface 
Approximation       Critical  radius 


boundary  conditions 
Per  cent  error 


Inward  flux 


10 


7.3976 
7.3201 
7.2805 
7.2826 
7.2784 
7.2793 
7.2778 
7.2783 
7.2775 
7.2778 

Mark's  vacuum-interface  boundary  conditions 


1.654 

1.895 

X  10 

.590 

1.251 

X  10 

.045 

5.433 

X  10 

.074 

6.973 

X  10 

.016 

3.221 

X  10 

.029 

4.370 

X  10 

.008 

2.216 

X  10 

.015 

3.056 

X  10 

.004 

1.S45 

X  10 

.008 

2.019 

X  10 

-3 


Critical 

Per  cent 

Weighted 

Per  cent 

Inward 

Approximation 

radius 

error 

average 

error 

f) 

ux 

•"l 

7.4979 

3.033 

8.644 

xlO-" 

7.3682 

1.250 

xlO-^ 

"2 

7.2359 

-.568 

2.140 

h 

7.2875 

.142 

3.076 

xlO-" 

7.2784 

.016 

xlO-^ 

•■4 

7.2663 

-.150 

1.204 

"•s 

7.2808 

.050 

1.791 

xlO-^ 

7.2776 

.005 

xlO-" 

^6 

7.2727 

-.062 

7.870 

'7 

7.2789 

.023 

1.241 

xlO-" 

7.2774 

.003 

xlO-" 

"8 

7.2749 

-.032 

5.650 

•"s 

7.2781 

.012 

8.147 

xlO-^ 

7.2773 

.001 

xlO-" 

P,« 

7.2759 

-.018 

4.167 

Table   II    (continued) 
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Infinite  black   reflector  boundary  conditions 

Critical          Per  cent          Weighted         Per  cent 
Approximation radius error average  error 


10 


7.5436 

3.661 

7.1960 

-1.116 

7.2962 

.261 

7.2562 

-.289 

7.2843 

.098 

7.2672 

-.137 

7.2810 

.052 

7.2714 

-.080 

7.2795 

.032 

7.2734 

-.052 

7.3716 


7.2791 


7.2773 


7.2772 


1.297 


.026 


.004 


.001 


Inward 
flux 

4 

.023 

X 

10 

2 

569 

X 

10 

1 

655 

X 

10 

1 

445 

X 

10 

6 

513 

X 

10' 

9 

950 

X 

lO' 

2 

949 

X 

10" 

7 

567 

X 

10" 

1 

760 

X 

10" 

6 

033 

X 

10" 

-3 


Approximation 


Variational   vacuum-interface  boundary  conditions 
Critical   radius  Per  cent  error 


7.3519 
7.2003 

7.2764 
7.2887 


1.026 
-1.057 
-.011 

.158 


Inward   flux 


2.374  X   10 


2.522   X   10 


7.010  X    10' 
1.643  X    10" 
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flux  at  tlie  vacuum  interface  over  all  angles  representing  neutrons  returning 
from  the  vacuum,  i.e.    f(.V-,u)'iu,   with  the  normalization  such  tiiat  tlie  total 
flux  at  tlie  center,  't'(O),  is  one  neutron  per  square  centimeter  per  second. 
Since  for  an  exact  solution  this  quantity  would  be  zero,  it  would  seem  that 
this  column  should  contain  numbers  which  would  be  indicative  of  the  accuracy 
of  the  critical  radius  computed  in  tliis  approximation.  However,  upon  close 
examination  it  is  readily  apparent  that  these  numbers  are  not  indicative  of 
the  accuracy  of  their  corresponding  critical  radii.  Even  for  a  particular 
boundary  condition  the  inward  flux  column  is  not  a  good  indication  of  the  accu- 
racy of  the  computed  critical  radii.  However,  for  a  particular  boundary  condi- 
tion the  nuraber  in  this  column  corresponding  to  odd  order  approximations 
decreases  for  increasing  orders  of  approximations.  Similarly  the  even  order 
numbers  decrease  with  increasing  orders  of  approximations  but  no  indication  of 
the  relative  accuracy  of  a  particular  computed  critical  radius  can  be  found  by 
comparing  these  numbers  for  the  even  and  odd  order  approximations.  The  basic 
problem  here  is  that  since  the  critical  radii  and  roots  for  each  of  the  approx- 
imations are  different  so  are  the  values  of  the  arguments  of  the  geometrical 
functions  and  hence  the  values  of  tlie  geometrical  functions  tliemselves.  Since 
this  column  has  little  indicative  ability,  it  will  be  omitted  from  the  remain- 
ing tables  of  results. 

Tlie  computed  results  with  Marshak  boundary  conditions  used  at  the  vacuum 
interface  of  a  bare  sphere  with  c  =  1.05  are  graphed  in  I'ig.  3.  This  type  of 
graph  where  the  computed  critical  radius  is  plotted  against  the  reciprocal  of 
the  order  of  approxiination  will  be  referred  to  as  a  convergence  graph. 
Although  the  graph  is  composed  of  discrete  points,  it  is  convenient  to  draw  in 
curves  connecting  the  points  representing  the  even  and  odd  order  approximations 
since  this  facilitates  visualization  of  the  rate  of  convergence.   It  is  readily 
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7.50 


7.40 


7.30 


ODD       ORDER       APPROXIMATIONS 
EVEN       ORDER        APPROXIMATIONS 


7.20  - 


0.0 


0.2  0.4  0.6  0.8 

RECIPROCAL       ORDER        OF        APPROXIMATION,       l/L 


1.0 
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apparent  from  Fig.  3  that  although  the  P,  result  is  more  accurate  than  the  Pj 
result,  the  odd  order  approximations  converge  to  the  exact  critical  radius 
more  rapidly  than  their  corresponding  even  order  approximations.   Both  the  even 
and  odd  order  approximations  tend  to  asymptotically  approach  the  exact  critical 
radius  from  above  along  separate  paths,  i.e.,  the  even  and  odd  order  approxi- 
mations are  not  counterconvergont.  Examining  tlie  Marshak  results  in  Table  II 
it  is  readily  apparent  that  as  the  order  of  approximation  increases  the  odd 
order  approximations  become  about  twice  as  accurate  as  their  successive  even 
order  approximations.  This  is  the  same  ratio  of  accuracy  of  odd  and  successive 
even  order  approximations  as  Davison's  formula  predicts  for  large  orders  of 
approximations,  in  spite  of  the  fact  that  when  using  Marshak's  vacuum- 
interface  boundary  conditions  the  even  and  odd  order  approximations  are  not 
counterconvergont  to  tlie  exact  critical  radius. 

Figure  4  is  a  convergence  graph  for  a  bare  spherical  reactor  in  which 
c  =  1.05  and  Mark's  boundary  conditions  are  used  at  the  vacuum  interface.  The 
weighted  average  points  in  Fig.  4  and  the  weighted  average  column  in  the 
appropriate  section  of  Table  11  are  found  by  applying  Eq.  (64)  with  the 
weighting  factors  given  in  Table  I.  The  expression  l''i  •'Y+i  Lw>  i"  which  L 
is  assumed  to  be  odd,  will  be  used  to  indicate  a  weighted  average  of  the 
critical  radii  from  the  P.  and  P.^,  approximations.   It  is  readily  apparent 
that  the  weighted  average  values  converge  to  the  exact  critical  radius  much 
more  rapidly  than  do  either  the  even  or  odd  order  approximations.  The  counter- 
convergence  trend  which  Eq.  (28)  predicts  for  Mark's  boundary  conditions  is 
also  shown  in  Fig.  4.  As  with  the  Marshak  boundary  conditions,  the  P.  result 
is  more  accurate  than  the  P  result,  but  for  higher  orders  of  approximations 
the  odd  order  results  are  better  than  tlie  corresponding  even  order  results. 
The  fact  that  the  weighted  average  critical  radii  converge  so  rapidly  to  the 
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exact  result  indicates  that  Davison's  formula  predicts  tlie  correct  trend  for 
this  case.  Since  tlie  |P.,P  |   is  not  as  accurate  as  might  be  expected,  the 
0(  "^  1*  }  t"^™  i"  ^'i-    (28)  is  not  negligible  for  tliese  low  orders  of 
approximation. 

Figure  5  is  the  convergence  graph  for  a  bare  spherical  reactor  with 
c  =  1.05  which  is  infinitely  reflected  with  a  black  medium.  As  with  Mark's 
vacuum-interface  boundary  conditions,  the  weighted  average  critical  radii 
converge  to  the  exact  critical  radius  much  more  rapidly  than  do  eitlier  the 
even  or  odd  order  approximations.   For  this  two  region  problem  the  counter- 
convergence  of  the  even  and  odd  order  approximations  is  again  readily  appar- 
ent. The  P.  result  is  once  again  more  accurate  than  the  P.  result  although 
for  higher  orders  of  approximations  the  odd  order  approximations  are  more 
accurate  than  the  corresponding  even  order  approximations.  The  iPiifilav 
for  tliis  case  is  not  as  accurate  as  would  be  expected  so  that  again  it  is 
apparent  that  the  higher  order  terms  in  Eq.  (28J  are  not  negligible. 

The  computed  results  for  a  bare  spherical  reactor  with  c  =  l.OS  using 
variational  boundary  conditions  at  the  vacuum  interface  are  shown  in  Fig.  6. 
From  the  limited  results  shown  it  is  not  possible  to  deduce  the  type  of  con- 
vergence which  the  variational  boundary  conditions  yield.   It  is  readily 
apparent  that  the  convergence  is  neither  counterconvergent  nor  asymptotic  as 
is  tlie  case  for  Mark's  and  Marsliak's  boundary  conditions  respectively.   In 
contrast  to  previous  cases  the  P.  and  V     approximations  with  variational 
boundary  conditions  yield  critical  radii  which  have  about  the  same  magnitude 
of  error  but  opposite  signs.   For  the  variational  boundary  conditions  the  P 
result  very  closely  approximates  the  exact  critical  radius  and  is  definitely 
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Figure  7  plots  tlie  angular  flux  at  the  vacuum  interface  as  a  function  of 
cose  for  a  bare  spherical  reactor  in  wliicli  c  =  1.05  resulting  from  a  F 
approximation  for  each  of  the  four  vacuum- interface  boundary  coniJitions  con- 
sidered. The  angular  fluxes  are  normalized  such  that  the  total  neutron  flux 
at  the  center,  ^CO),  is  one  neutron  per  square  centimeter  per  second,   for  an 
exact  solution  Eq.  (29)  shows  tliat  in  the  range  -1  <^  u  f_  0  the  angular  flux 
would  be  zero.   I'igure  7  shows  that  tlie  variational  boundary  conditions  best 
approximate  this  condition.   It  is  also  apparent  from  Fig.  7  and  Table  I  that 
the  closeness  with  which  the  various  boundary  conditions  approximate  the 
exact  condition  the  closer  the  critical  radius  is  to  the  exact  critical  radius. 
However,  the  relative  accuracy  with  which  a  particular  boundary  condition 
approximates  the  exact  condition  is  not  directly  proportional  to  the  relative 
accuracy  of  the  computed  critical  radius.  The  type  of  graph  shown  in  Fig.  7 
could  be  employed  as  a  rough  guide  to  find  out  which  critical  radius  deter- 
mined from  a  set  of  boundary  conditions  is  the  most  accurate  for  a  particular 
problem. 

Tables  III  and  IV  list  the  computed  results  for  a  bare  spherical  reactor 
in  whicli  the  values  of  c  are  1.02  and  1.4U  respectively.  As  the  value  of  c 
deviates  from  unity  by  increasing  amounts,  it  seems  natural  tliat  more  and  more 
terms  in  tlie  Legendre  polynomial  expansion  of  tlie  angular  flux  will  be  needed 
to  satisfy  the  exact  boundary  condition  to  a  preset  degree  of  accuracy.  Thus 
it  will  be  expected  that  the  rate  of  convergence  of  the  spherical  liarmonics 
approximations  will  be  proportional  to  the  relative  departure  of  the  value  of 
c  from  unity.  This  trend  is  readily  apparent  in  the  per  cent  error  columns  of 
Tables  II,  III,  and  IV. 

Some  generalizations  can  be  made  as  to  the  relative  accuracy  obtained 
with  various  boundary  conditions  by  examining  Tables  II,  III  and  IV.  The 
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Table  III.  Calculated  results  for  a  bare  spherical  reactor  in  which  c 
is  1.02.  The  exact  critical  radius  is  12.0270  mean  free  paths. 


Marshak  vacuum-interface  boundary  conditions 
Approximation Critical  radius Per  cent  error 


12.1269 
12.0729 
12.0312 
12.0338 
12.0289 
12.0301 
12.0282 
12.0288 
12.0279 
12.0283 


0.831 
.382 
.035 
.057 
.016 
.026 
.010 
.015 
.007 
.011 


Mark's  vacuum-interface  boundary  conditions 


Approximation 

Critical 
radius 

Per  cent 
error 

Weighted 
averaqe 

Per  cent 
error 

12.2239 

11.9787 

1.637 
-.402 

12.1025 

0.628 

"3 
•■4 

12.0395 
12.0138 

.104 
-.110 

12.0285 

.012 

12.0319 
12.0214 

.041 
-.047 

12.0277 

.006 

12.0297 
12.0242 

.022 
-.023 

12.0276 

.005 
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Table  III  (continued) 


Infinite  black  reflector  boundary  conditions 

Critical    Per  cent    Weighted     Per  cent 
Approximation radius error average error 


12.2520 

1.871 

12.1060 

11.9571 

-.581 

12.0448 

.148 

12.0289 

12.0076 

-.161 

12.0343 

.061 

12.0278 

12.0180 

-.075 

12.0312 

.035 

12.0276 

12.0219 

-.042 

12.0298 

.023 

12.0276 

12.0238 

-.027 

0.657 


.016 


.007 


.005 


.005 


Variational  vacuum-interface  boundary  conditions 
Approximation Critical  radius Per  cent  error 


12.0827 
11.9843 
12.0262 
12.0336 


0.463 

-.355 

-.007 

.055 


Table  IV.  Calculated  results  for  a  bare  spherical  reactor  in  which  c 
is  1.40.  The  exact  critical  radius  is  1.9854  mean  free  paths. 


Marshak  vacuum-interface  boundary  conditions 
Approximation Critical  radius Per  cent  error 


2.1225 
1.9760 
1.9869 
1.9873 
1.9860 
1.9861 
1.9856 
1.9856 
1.9855 
1.985S 


6.895 
-.473 
.076 
.096 
.030 
.035 
.010 
.010 
.005 
.005 


Mark's  vacuum-interface  boundary  conditions 


Critical 

Per  cent 

Weighted 

Per 

cent 

Approximation 

radius 

error 

average 

error 

h 

2.2224 

11.937 

2.0920 

5 

.369 

h 

1.9590 

-1.330 

h 

1.9898 

.222 

1.9865 

.055 

'"4 

1.9821 

-.166 

■"5 

1.9867 

.065 

1.9857 

.015 

"6 

1.9843 

-.055 

''7 

1.9859 

.025 

1.9855 

.005 

'■8 

1.9849 

-.025 

"9 

1.9856 

.010 

1.9854 

0 

P,n 

1.9851 

-.015 

Table  IV  (continued) 


Infinite  black  reflector  boundary  conditions 

Critical    Per  cent    Weighted     Per  cent 
Approximation radius error average      error 


10 


2.3531 

18.520 

1.7629 

-11.207 

2.0394 

2.720 

1.9430 

-2.136 

1.9994 

.705 

1.9692 

-.816 

1.9914 

.302 

1.9766 

-.443 

1.9889 

.176 

1.9798 

-.282 

2.0610 


1.9982 


1.9855 


3.808 


.645 


.096 


.005 


Variational  vacuum-interface  boundary  conditions 
Approximation Critical  radius Per  cent  error 


1 


2.0779 
1.7481 
1.9858 
2.0605 


4.659 

-11.952 

.020 

3.783 
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spherical  liarmonics  approximations  are  seen  to  converge  iiiucli  more  rapidly 
for  Mark's  boundary  conditions  than  for  the  infinite  black  reflector  boundary 
conditions.   I'or  slab  geometry  tliese  two  sets  of  boundary  conditions  are 
equivalent,  but  in  spherical  geometry  they  are  seen  to  be  different  although 
tliey  botli  counterconverge  to  the  exact  critical  radius  as  expected.  The 
difference  in  spherical  geometry  is  undoubtedly  due  to  the  fact  that  the  infi- 
nite black  reflector  boundary  conditions  take  the  curvature  of  the  outer  edge 
of  the  spherical  reactor  into  account  in  a  much  different  manner  than  do 
Mark's  boundary  conditions.  Lixcept  for  the  infinite  black  reflector  boundary 
conditions,  the  boundary  conditions  used  in  this  work  are  actually  derived  for 
slab  geometry  cases.  The  fact  tiiat  these  slab  geometry  boundary  conditiojis 
yeild  such  good  results  indicates  that  tliese  boundary  conditions  are  also 
applicable  in  spherical  geometry  where  the  curvature  of  the  vacuum  interface 
is  not  too  great. 

It  was  previously  stated  that  the  variational  boundary  conditions  would 
be  expected  to  give  tlie  most  accurate  results  since  they  are  developed  directly 
from  the  matliematics.   It  is  readily  apparent  in  iables  11,  III  and  IV  tliat, 
at  least  for  the  1'  and  P  approximations,  tiu-  variational  boundary  conditions 
give  the  most  accurate  critical  radii  for  those  orders  of  approximations.   No 
particular  set  of  boundary  conditions  yields  a  best  critical  radius  for  tlie  I' 
approximation.   As  expected,  Marshak's  boundary  conditions  are  seen  to  be  loss 
accurate  tlian  the  variational  boundary  conditions.   In  addition,  for  a  given 
order  of  approximation,  Mark's  or  the  infinite  black  reflector  boundary  condi- 
tions are  always  poorer  tiian  Marshak's  boundary  conditions.   However,  the 
weighted  average  critical  radii  from  Mark's  boundary  conditions  arc  seen  to 
converge  more  rapidly  to  tlio  exact  critical  radius  tlian  tlie  odd  order  approxi- 
mations using  Marsliak's  boundary  conditions,   hxcept  for  tlie  1 1'  1'  I    the 

I    1'  2'av' 
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tions  are  about  as  accurate  as  the  P    Marshak  boundary  condition  results, 

e.g.,  the  |P,,P,|   using  Mark's  boundary  conditions  is  about  as  accurate  as 
"   '  '  3  4 ' av     "  ' 

the  Pg  Marshak  boundary  condition  result.  The  critical  radii  resulting  from 


about  as  accurate  as  the  iPr,'',  I   from  Mark's  or  the  infinite  black  reflector 
'  5 '  0 '  av 


conditions  applied  at  the  vacuum  interface  yields  the  most  accurate  estimate 
of  the  critical  radius. 
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3.3  kesults  of  Numerical   Calculations 
for  Ueflectui)  Splierical   Koactors 

In  order  to  examine  tlie  rate  of  convergence  of  the  spherical  liarmonics 
approximations  for  a  reflected  spherical  reactor  the  cases  shown  in  Table  V 
will  be  considered. 

Table  V 
Reflected   Ueactor  Cases   to  be  Considered 
Case  c  1/1  RgCmfp) 

1  O.yj  1.0  1.0 

2  .99  l.U  10.0 

3  .50  1.0  2.0 

4  .99  2.5  1.0 

5  .99  2.5  10.0 

6  .50  2.5  2.0 

The  combinations  of  Tr  and  c^  for  the  six  cases  are  such  that  for  cases  1  and 

o 

4  the   reflector  is   about   0.2  diffusion   lengths   in  tliickness  whereas   for  cases 
2,   3,   5,   and  b   tlie   reflector  is   about   2.0  diffusion   lencths   in  thickness.     As 
previously  noted,   a  reflector  wliicli   is   1.5   to  2.0  diffusion   lengths   in  thick- 
ness  is  essentially   infinite.      iVhen  a  reflector  is  of  sucli   a  thickness,    it  will 
be   referred  to  in  this  work  as   a  "thick"  reflector.  .  In  contrast  when  tlie 
thickness   of  the    reflector   is  much    less   than   2.0  diffusion   lengtlis,    it   will   be 
termed  a  "thin"   reflector.      Thus    cases    1    and   4   represent   "tliin"   reflectors 
whereas   cases  2,    3,   5,    and  (j   represent   "tliick"    reflectors.      Tlie   value  of 
c  =    .99   represents   a   reflector  which    is   a   good  moderator  whereas   tlie   value  of 
c  =   .50   represents   a   reflector  wliich   is   a  breeder  blanket   or  other  poor 


moderator.     The  case  of  a  small  value  of  c   for  a  thiri  reflector   is  not   con- 
sidered since   for  sucli  a  case  tlie   tliickness  of  the   reflector  would  1)0  about 
0.2  mean   free  paths  or  less   and  hence  would  be  extremely  thin  physically.      In 
cases   1,   2,   and  3  the  total  macroscopic   cross   section,  Z,    is  continuous  across 
tlie  core-reflector  interface  whereas   in  cases   4,   5,   and  6  the  total  macro- 
scopic cross   section   is  discontinuous   across  the  boundary. 

As   in  the  bare   spherical   reactor  cases,    l.OS  will  be  used  as  the  value  of 
c,  the  mean  number  of  secondaries,    in   a  basic  reflected  reactor  core  to  which 
the   results   for  tlie  values  of  c  =   1.02   and   1.40  will  be  compared.     Whenever  a 
counterconvergent  trend  is   found  in  the  computed  results,   Uq.    (64)   will  be  used 
to  find  weighted  average   critical   radii   and  they  will  be   listed  in  a  weiglited 
average  column.     Table  VI    lists   the  computed   results   for  case   1  with  c  =   1.05 
in  the  core  region.      It   is  readily  apparent   from  the  data  that   the  nature  of 
tlie  convergence   for  each  of  the  vacuum-interface  boundary  conditions  is  the 
same   as   the   convergence  with   the   boundary  condition   applied  to   a  bare   spherical 
reactor.     Thus   for  a  thin   reflector  and  continuity  or  near  continuity  of  the 
total  macroscopic  cross   section,   and  the   mean  number  of  secondaries  across 
tl»e  coro-ref lector  boundary,   the  vacuum-interface  boundary  conditions   dominate 
the   convergence  pattern.      Tables   VI 1    and   VIll    list   the   computed  results    for 
cases  2   and  3   respectively  with  c  =   1.05   in  the   core   region.     The   convergence 
pattern   in  these  two  cases   is   seen  to  reflect  the  counterconvergent  trend 
expected  of  an   interfacial   boundary  or  Mark's  vacuum- interface  boundary  condi- 
tions.    Thus   for  a  thick   reflector  and  continuity  of  the  total   macroscopic 
cross   section   across   the   core-reflector  interface  the    interfacial   boundary 
dominates  the   convergence  jiattern.      The  computed   results    for  a  reflected 
reactor  whose  constants   are  tliose  of  case   4   and  whose   core  has   a  value  of  c 
equal  to   1.05   are   listed   in  Table   IX.      li   is   readily  apparent   from  the   results 
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Table  VI.      Calculated  results  I'or  case   1    (reflected  reactor) 
in  wliicli  c  is   1.05. 


Approximation 


Marsliak's  B.C. 


Critical   radius 


Mark's  B.C. 


Variational  B.C. 


Critical    Weighted 
radius average Critical  radius 


6.6192 

6.6901 

6.S114 

6.4600 

6.4873 

6.4910 

6.4880 

6.4788 

6.4862 

6.4876 

6.4863 

6.4824 

6.5762 


6.4858 


6.5875 
6.4274 
6.4852 
5.4943 


Table  VII.      Calculated  results   lor  case  2   (reflected  reactor) 
in   wliicli  c   is    1.05 


Marsliak's  B.C. 


Mark's   B.C. 


Variational   B.C. 


Critical       Weiglited       Critical       Weiglited       Critical       Weighted 
Approximation       radius average         radius average         average         average 


5.1896 
5.0625 
5.0779 
5.0739 
5.0762 
5.0749 


5.1913 

5.0609 
5.0780 
5 . 0736 
5.0763 
5.0748 


5.1268 


5.0761 


5.0757 


5.1889 
5.0611 
5.0778 
5.0740 


5.1256 


5.0762 
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Table  VIII.     Calculated  results  for  case  3   (reflected  reactor) 
in  which  c  is   l.OS. 


Marsliak's   B.C. 


Mark's   B.C. 


Variational   ii.C. 


Critical       Weighted       Critical       Weighted       Critical       Weighted 
Approximation       radius average         radius average         radius average 


7.3215 
7.0284 
7.0981 
7.0731 
7.0903 
7.0803 


7.1764 


7.0874 


7.0863 


7.3223 
7.0283 
7.0983 
7.0733 
7.0903 
7.0803 


7.1768 


7.0876 


7.0863 


7.3211 
7.0271 
7.0981 
7.0734 


7. 1756 


7.0875 


Table  IX.  Calculated  results  for  case  4  (reflected  reactor) 
in  which  c  is  1.05, 


Marsliak's  B.C. 


Mark's  B.C. 


Variational  B.C. 


Critical   Weighted   Critical   Weighted   Critical   Weighted 
Approximation   radius average    radius average    radius average 


6.9287 
6.6506 
6.6967 
6.6829 
6.6942 
6.6877 


6.7911 


6.6908 


6.9788 
6.6370 
6.6981 
6.6813 
6 . 6946 
6.6864 


6.6909 


6.9062 
6.5510 
6.6965 
6.7025 


6.7304 


6.6991 
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for  tliis   tliin   reflector  case  that   the  couiiterconvergunt  trend   imposed  by  tlie 
core-rof lector  interface   is   tlie  predominant   factor  in  the   convergence  pattern 
although  the  vacuum-interface  boundary  conditions  liave  considerably  more 
effect   on  the  convergence  pattern  than  was  observed  in  tlie  two  previous  thick 
reflector  cases,      llms,    for  a  thin  reflector  and  discontinuity  of  the  total 
macroscopic  cross  section   at   the  coro-reflector  interface  the   intorfacial 
boundary  conditions  dominate  the  convergence  pattern  although  the  vacuum- 
interface  boundary  conditions  may  play  a   significant   role   in  the   relative 
speed  of  convergence.     Tables   X  and   XI   list   the   computed  results   for  cases  5 
and  6  for  a  reflected  reactor  wliich  lias   c  =   l.OS   in  the  core   region.     Here  as 
with  the  other  thick   reflector  cases  the  core-reflector  interface  dominates 
the  convergence  picture.      For  thick   reflectors  the   interfacial  boundary  domi- 
nates the  convergence  pattern  so  strongly  tliat   the  vacuum-interface  boundary 
conditions  have   little   or  jio  effect  on  the   convergence.     Thus,   as  would  be 
expected,   the   convergence  trend   for  a  reactor  reflected  with   a  thick   reflector 
is   always  that   of  counterconvergence.      Tables   XII    and  XIII    list    the   computed 
results   for  case  5  with  the  values   of  c  in  the  core  equal  to  1.02  and   1.40 
respectively.     As  with  the  previous  thick   reflector  cases,   the   interfacial 
boundary  dominates   the  convergence  pattern.      It   is   readily  apparent   that   when 
the  core   region  has   c  =   1.40  the  vacuum-interface  boundary  conditions  have 
little  or  no  effect   on   tlie   convergence.      In   general   the  greater  the  disconti- 
nuity in  the  value  of  c  at   the   core-reflector  interface  the  greater  the 
dominance  of  the   interfacial   boundary  in  the   convergence   pattern. 

In  all   of  the  computed  results   for  the   reflected  reactor  cases  tiie 
weighted  average   critical   radii,   whenever  noted   in  the   tables,   converge  more 
rapidly  than  the  critical   radii   computed  directly   from  the   spherical  harmonics 


Table   X.      Calculated   results   for  case  5    (reflected  reactor) 
ill  which  c   is   1.05. 


Marshak's  B.C.  Mark's   B.C.  Variational   B.C. 

Critical       Weiglited  Critical       Weighted       Critical       Weighted 

Approximation       radius average  radius average         radius average 

P  6.1999  6.2008  6.1995 

'  0.0398  6.0399  6.0392 

P  5.8765  5.8757  5.8757 

P  5.9697  5.9698  5.9697 

^  5.9577  5.9577  5.9577 

P,  5.9416  5.9414  5.9416 

4 

P  5.9580  5.9581 

5.9544  5.9544 

P^  5.9489  5.9488 


Table  XI.     Calculated  results   for  case  6   (reflected  reactor) 
ill  which  c  is   1.05. 


Marshak's   B.C.  Mark's   B.C.  Variational   B.C. 

Critical       Weiglited       Critical       Weighted       Critical       Weighted 
A|)proxiiiiation       radius average         radiiiii  average  radius average 


''l 

7.4274 

'2 

6.9812 

''3 

7.1508 

''4 

7.0900 

''5 

7.1306 

''6 

7.1069 

7.1248 


7.4281  7.4271 

7.2071 

6.9816  6.9792 

7.1511  7.1507 

7.1252 

7.0905  7.0900 


7.2054 


7.1307 

7.1211  7.1212 

7.1069 


* 
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Table   XII.     CalculateJ  results   for  case  5    (reflected  reactor) 
in  wliich  c   is    1.02. 


Marsliak's   b.C.  Mark's   B.C.  Variational   li.C. 

Critical       Wciglited       Critical       Weighted       Critical       WeiRhted 
Approximation        radius ave  ra  .i;o radius averaj^e radius average 

10.2109 

10.1068  10.1057 

9.9984 


■"l 

10.2116 

10.2133 

10 

1068 

h 

9.9998 

9.9982 

"3 

10.UJS9 

10.0360 

10 

0314 

"4 

10.0254 

10.0251 

''5 

10.UU8 

10.0319 

10 

U303 

1' 

10.0281 

10.0279 

10.0313 


10.0303 


10.0358 
10.0254 


Table   Xlll.      Calculateil   results    for   case   5    (reflected   reactor) 
in   wliich   c    is    1.40. 


'larshak's   B.C.  Mark's   U.C.  Variatiojial    B.C. 

Critical       Weialited       Critical       Weiglited       Critical       Weighted 
Approximation        radius average radius average radius average 

2.0659  2.0657 

1.5012  1.5013  1.5012 

.9252  .9252 


''l 

2.0657 

''2 

.9252 

"3 

1.8000 

% 

1.5200 

h 

1.7265 

't 

1.6187 

1.6803 


1.8000  1.8000 

1.6803 
1.5200  1.5200 


1.7265 
1.6187 


1.6834 
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approximations  except  that  the  jP-.F, |,   is  usually  not  any  better  than  the 
Pj  result.   These  observations  are  in  agreement  with  tlie  bare  splierical 
reactor  cases  for  Mark's  boundary  conditions  or  the  infinite  black  reflector 
boundary  conditions. 

Ill  the  cases  considered  here  it  is  found  tliat  the  spatial  dependence  of 
tlie  total  neutron  flux  is  similar  to  that  of  Fig.  2.   Figure  8  is  a  plot  of 
tlie  angular  flux  at  the  core-reflector  interface  as  a  function  of  cos9  as 
viewed  from  both  sides  of  the  boundary  for  a  V     approximation.   It  is  apparent 
that  there  is  a  discontinuity  in  the  angular  flux  across  the  boundary  near 
cose  =  0  but  that  tlie  dependence  of  the  angular  flux  on  the  radial  variable 
is  more  nearly  continuous  as  cose  approaches  the  value  of  '1.  Tliis  means  that 
the  normal  component  of  tlie  angular  flux  is  nearly  continuous  at  the  inter- 
facial  boundary  but  that  the  tangential  component  is  discontinuous.  This 
condition  arises  from  the  fact  that  in  an  even  order  approximation  the  expres- 
sion f.W   -   P,(0)(2)t+l):f.,(K)  is  matched  at  the  interfacial  boundary.  This 
expression  can  be  shown  (3)  to  be  equivalent  to  removing  the  dependence  of 
/(r,w)  O"  r  for  small  values  of  p(cose);  thus  allowing  a  discontinuity  in 
/(K.m)  for  small  values  of  u.  The  reason  wh/  fiT,\i)    is  not  continuous  across 
the  interface  for  values  of  jj  approaching  unity  is  that  a  finite  instead  of 
infinite  order  of  approximation  is  being  utilized.   A  similar  plot  for  an  odd 
order  approximation  has  the  two  curves  in  Fig.  8  coinciding  since  continuity 
of  the  angular  flux  at  the  interfacial  boundary  in  an  odd  order  approximation 
is  required  by  Eq.  (25). 

From  the  preceding  discussion  of  the  computed  results  for  reflected 
reactors  three  factors  appear  to  be  important  in  the  convergence  pattern. 
First,  the  relative  thickness  of  the  reflector  in  diffusion  lengths  is  impor- 
tant.  For  thick  reflectors  the  core-reflector  interface  dominates  the 
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FROM      CORE      SIDE 
FROM       REFLECTOR       SIDE 


-1.0 


-0,5 


0.0 
COS  e 


0.5 


1.0 


Figure     8.     ANGULAR      FLUX      AS      A      FUNCTION     OF      COS  6 
AT     THE      CORE  -  REFLECTOR      INTERFACE     FOR 
A     P^      APPROXIMATION. 
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convergence  pattern  and  produces  a  counterconvergent  trend.     The  second  and 
third  important  factors  in  the  convergence  pattern  are  the  relative  size  of 
the  discontinuities   in  the  values  of  the  mean  number  of  secondaries,   and  the 
total  macroscopic  cross  section  at  the  core-reflector  interface.     When  a  sig- 
nificant discontinuity  in  either  c  or  Z  exists  at  the  core-reflector  interface 
the  convergence  pattern  is  one  of  counterconvergence,  thus  demonstrating  that 
the   interfacial  boundary  is  the  dominant   factor  in  tlie  convergence  pattern 
for  such  cases  even   for  a  relatively  thin  reflector. 
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3.4  Conclusions 

The  A     constant  was  assumed  to  be  zero  in  computing  the   critical   radius 
of  a  bare  splierical   reactor  when   all  of  the  vacuum-interface  boundary  condi- 
tions except  the  infinite  black   reflector  boundary  conditions  were  used.     The 
consequences  of  this   assumption  will  now  be  considered.      Both  Mark's  vacuum- 
interface  boundary  conditions  and  the  infinite  black  reflector  boundary  condi- 
tions brought  about  convergence  trends  which  counterconverged  in  a  like 
manner.     Since  the   infinite  black   reflector  boundary  conditions  take  A     into 
account  whereas  Mark's  boundary  conditions  do  not,   the  assumption  that  A     was 
zero  at  a  vacuum  interface  seems  entirely  justified. 


of  the  exact   critical   radius.      For  higher  orders  of  approximations  the  odd 
order  approximations  were   invariably  found  to  be  superior  to  their  correspond- 
ing even  order  approximations.      Thus  the  popular  belief  that  the  even  order 
approximations  are  inferior  to  the  odd  order  approximations  is  not  entirely 
justified. 

When  Davison's   formula,   which  predicts  the  accuracy  of  the   L       order 
approximation,  was   rearranged  to  yield  weighting  factors,   it  was  seen  to 
predict  very  accurate  weighted  average  critical  radii   for  the   P,   and  higher 
order  approximations.      Althougli  Davison's   formula  was  developed   for  slab 
geometry  it  was   seen  to  be  quite  valid  in   spherical   geometry  as  well.      For  the 
P.   and  P     approximations  Davison's   formula,   as  used  here,   was  definitely 
incorrect   indicating  that   for  these   low  order  approximations  the  higher  order 
terras   in  Eq.    (28)   cannot  be  ignored. 
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tional  boundary  conditions  yielded  a  critical   radius  which  was  about   as  accu- 
rate as  a   IP^.P. I        usinc  Mark's  or  the   infinite  black  reflector  boundary 
'    5*   6'av 

conditions,   or  the   V     result  with  Marshak's  boundary  conditions.     The 
|P,,P,^j|        for  Mark's  or  the  infinite  black  reflector  boundary  conditions 
yielded  critical  radii  which  were  of  the  same  order  of  accuracy  as  a  more 
involved  P         approximation  using  Marshak's  boundary  conditions.      This   last 
observation  was   found  to  be  invalid  for  the   P     and  P.   approximations.      For  the 


boundary  condition  consistently  yielded  the  best  estimate  of  the  exact   critical 
radius. 

For  reflected  reactor  cases  a  counterconvergent  pattern   in  the  conver- 
gence,  due  to  the  core-reflector  interface,   was  predominant  whenever;   a.)   the 
reflector  was  "thick";   b.)   there  was  a   large  discontinuity  in  the  value  of  the 
total  macroscopic  cross   section,   £,   across  the   core-reflector  interface;   and/or 
c.)   there  was   a  large  discontinuity  in  the  mean  number  of  secondaries  per 
interaction,   c,   across  the   core-reflector  interface.      For  all  but  the   P     and 
P     approximations  the  weighted  average   critical   radii   for  reflected  reactors, 
in  which  the  core-reflector  interface   dominated  the  convergence  pattern, 
appeared  to  be  closer  to  the  exact  critical   radius  than  the  results   from  either 
the  even  or  odd  order  approximations. 
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4.0  SUGGESTIONS   FOR  FURTIini!  STUDY 

In  this   study  only  a  few  representative  values  of  c  are  chosen   for 

investigation  of  the  nature  of  the  convergence  of  the   spherical  harmonics 

approximations.      For  large  order  of  approximations   (three  or  greater)   when 

Davison's   formula  is  used  to  find  the  weighting  factors  in  Fq.    (65)   the 

weighted  average  critical  radii  appear  to  be  very  accurate.     Unfortunately, 

the  saine  cannot  be  said  for  tlie    IF,  ,1'   I      .      A  study  should  be  made  of  the 

'  1  2 ' av        ' 

higher  order  terras  in  Davison's  formula  so  that  it  could  be  improved  for  low 
order  approximations.   Alternatively  a  study  could  be  made  to  determine  an 
empirical  formula  which  would  predict  the  values  of  the  weighting  factors 
over  a  wide  range  of  values  of  c.   If  either  of  these  studies  were  successful 
to  any  substantial  degree,  only  tlie  relatively  simple  P.  and  P  approximations 
would  be  needed  to  obtain  a  trustworthy  critical  radius. 

The  work  done  here  can  be  extended  to  the  simpler  case  of  slab  geometry 
by  noting  that  the  only  difference  between  the  slab  and  spherical  geometries 
is  the  difference  between  the  Q.(x)  and  e"  geometrical  functions.   This  exten- 
sion would  only  involve  some  small  changes  in  the  computer  program  described 
in  Appendix  C,   Also  this  work  should  be  extended  to  multigroup  approxima- 
tions.  It  is  possible  that  many  parts  of  the  program  described  in  Appendix  C 
could  be  used  in  such  a  study. 

Since  the  variational  boundary  conditions  were' found  to  be  the  most  accu- 
rate of  all  the  boundary  conditions  considered,  it  is  logical  to  suggest  that 
they  be  extended  to  approximations  higher  than  the  P  and  P  approximations 
used  here  if  accuracies  greater  than  those  resulting  from  the  P,  approximation 
are  required.   The  odd  order  variational  boundary  conditions  have  been  found 
to  be  more  accurate  tlian  the  corresponding  even  ordor  variational  boundary 
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conditions.  Therefore  the  most  fruitful  results  will  probably  be  obtained  by 
using  the  odd  order  variational  boundary  conditions. 
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APPENDIX  A 


Tlie  Q,(x)    function  is  defined  by  Weinberg  and  Wigner   (21)   as 


where  ''.^.i/Tt*^  i^  *  modified  Bessel  function  of  the  second  kind.  The  Q,(x) 
functions  are  2/ii  times  the  modified  spherical  Bessel  functions  of  the  third 
kind.     Watson   (19)   defines   tlie   K     ,/,(x)   functions  as 

"4+1/2 '"' °  "27®    ^  r-  f*"^^ 

''*^'^  ^^  j=0  j!(«-j)!(2x)5 


explicitly  as 

X  £ 


(t*3)l 


Q,(x)    -  I    I       '-"    ''''>•■    ,  Q,(-x)    --'-I  ''-'"        .    .        (A-3) 

*  j=0  j!(£-j)l(2x)^  "  "     j=0  jl(4-j)!(2x)J 


Q„(X)     =    2.      ;    Q    (X)     =    2.    (1    -    i)     ;    Q    (x)     =    £    (1    -    1   ♦   i^) 


'^tt")    =   V2(''i    -^VlW    •  fA-") 


and 


(a7^^)*^(='^  =  ViW  (A-5) 

follow  immediately  from  the  properties  of  Bessel   functions. 
In  this  work  a  new  function,   C   (x)   was  defined  as 

(A-6) 
The  C   (x)    functions  arc  tlie  modified  spherical   Bessel   functions  of  the   first 
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kind.  Mark  (10,11)  considers  these  functions  and  shows  that  for  small  x. 


{1    +   0(x^))    .  (A-7) 


'^i^"'    '   l-3-5---(2U'IT 
At  this  point   it   is  necessary  to   investigate  the   finiteness  of  these 
functions  at  x  =  0.     Since 

Um  C   (x)   =   am  YTY.s^.. (ii.i)    ^   *  Of''^) '   '  '^iC''   =   (J'.   I  >  0.        '■^'^^ 
x-K)  x->0 

the  C  (x)  functions  are  thus  bounded  for  a  zero  argument  and  so  the  total  neu- 
tron flux  is  bounded  at  r  =  0  as  is  required  in  section  2.3. 

The  C  (x)  functions  may  be  written  as  a  combination  of  sinh(x)  and 
coshCx)  terms.   In  this  form  the  C  (x)  functions  are,  for  I   even 

C  (X)  =  ""''t'^)  I  (J^*^)'     .  g°s''(x)  ^     (3*^'     ,     (A-9) 

"^        *   j=0  jl(i-j)l(2x)^      "   j=l  j!(i-j)l(2x)^ 
j  even  j  odd 

and  for  ii.  odd 

C  (x)  =  ''°^^'M    I  (J*M'     .  sinh(x)  ^     (J^lt)'     .    (A-10) 

"        "   j=0  j!(£-j)!(2x)'      *   j  =  l  jl(ll-j)!(2x)' 
j  even  j  odd 

The  arguments  of  the  C  (x)  functions  which  are  encountered  in  this  work  are 

either  entirely  real  or  entirely  imaginary.   Equations  (A-3) ,  (A-9),  and 

(A-10)  are  appropriate  for  real  arguments.   The  imaginary  arguments  arise  only 

in  the  central  core  region  in  which  tlie  C  (x)  functions  are  used.   In  this 

region  one  of  the  roots  X,  is  imaginary.   The  argument  for  the  C  (x)  functions 

is  tr/X.,   so  the  form  of  the  imaginary  argument  is  x/i  or  -ix.   For  this  type 

of  imaginary  argument  the  C  (x)  functions  become,  for  It  even 


^        1     Mark  actually  defines  a  function  H.  (x)  wliich  is  the  same  as  Q|^(x)  in  this 
development  and  considers  the  combination  H.  (x)  *   (-1)  "|j(-x)  which  differs 
from  the  C.  (x)  used  here  only  by  a  factor  ot  two. 

■I 


i.^  '         If     t,  1 


and  for  £■  odd 


3=1   j!(£-j)!(2x)J 
j  odd 
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cos(x)  ^  (■i)C3-l)/2(^^j),  ^      ^^^^^^ 

"  j=l   j!()l-j)!(2x)' 
j  odd 


cos(x)  ^  (-t)J^^(Uj)!|     (A- 12) 
"  J=U  j!(t...i)!(2x)^ 


]  even 


2 
and  an  imaginary  root.   Since  -i  =  +1,  the  product  of  ti  (iX.)  and  C  (Sr/iX,  ) 


for  odd  1  will  be  positive  and  real. 
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APPENDIX  B 

Variational  Boundary  Conditions 

Pomraning  and  Clark  (17)  have  used  the  variational  calculus  to  develop  a 
set  of  consistent  vacuum-interface  boundary  conditions.  The  fundamental  equa- 
tion of  the  variational  boundary  conditions  at  a  right  hand  vacuum-interface 
boundary  is 


r. 


/(R.M)6/(«,-M)pdM  -  0  .  •  (B-1) 


In  order  to  arrive  at  a  set  of  coefficients  b  .  which  can  be  put  into  the  form 

''J 

L 

I  fAm,,    =0     J  =  l,2....,[(Ln)/2]  ,  (B-2) 

the  angular  flux  will  be  expanded  with  the  same  Legendre  polynomial  expansion 
that  Pomraning  and  Clark  used,  namely: 

/(R.m)  =  I  (^)»,(«)P,(M)  .  '      (B-3) 

1=0 

instead  of  the  expansion 

L 

fvi.v)  =  I  f,mi',M  (B-4) 

i'O   "    " 
which  is  used  throughout  most  of  this  work.   From  Eqs.  (B-3)  and  (B-4)  it  is 
recognized  that 

h^''^  -  rrn:  ^^^^  ■  ^^-'^ 

Pomraiiing  (16)  has  noted  that  in  order  to  solve  Eq.  (B-1)  there  must  be 


combinations  can  bo  put  into  the  general  form 

*2l(UlJ/2)-,n('*^  *  J„   \Jl^''^    -   "   f°^  1  im^  [(Ul)/2]  .  (B-6) 
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For  a  F'  approximation  this  relationship  reduces  to 

Substituting  Eq.  (B-3)  into  Eq.  (B-I)  and  using  Eq.  (B-7)  in  the  result,  Eq. 


i^-B-^oi  1  *of''^**oW  =  °  •  ^'-«^ 


Eq.    (B-8)   can  be   an  equality  if  and  only  if  the  coefficient  of  ♦yCRjS'tigCR)   is 
set  to  zero,   in  which  case 

'^01  =  2/a  , 

Uyj   -   *    /2/1   .  (B-9) 

Physical  considerations  demand  that  the  negative  sign  of  Eq.  (B-9)  be  used. 
Considering  Eqs.  (B-5),  (B-7),  and  (B-9)  it  is  readily  apparent  that  the 
coefficients  b  .  in  Eq.  (U-2),  for  a  P  approximation,  are  given  by 

X.J  1 

^01  =  1  .    "11  =  -  •^- 


tional  equation  comes  directly  from  the   spherical  harmonics  form  of  the  Boltz- 
inann  transport  equation.     Using  Eq.    (B-3)   ifi:itead  of  Eq.    (B-4)   to  expand  the 
angular  flux,   tlie  spherical  harmonics   form  of  the  Boltzmann  neutron  transport 
equation   in   slab  geometry,   subjected  to  tlie  restrictions  used  throughout  this 
work,   becomes 

2OT  37  n*l^^^    *   2OT  4  ♦.-l^''^    *   I(l-C6^^)*^(x)    =   0   .  (B-10) 

For  a  Pj   approximation  these  equations  become 

jjL4,j(x)    *   I(l-c)*y(x)    =   0    ,  (B-11) 
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l^'fiCx)  *   i:*2(x)  =  0  .  (B-13) 

From  Eqs.   (B-1)  and  (li-13)  the  elimination  of  the  ((),(x)  dependence  shows  that 

*2M   =  ^   (1-c)*q(x)  .  (B-14) 

Tliis   is  the  required  extra  relationship.     Now,   substitutinR  Eq.    (B-3)   into 
liq.    (U-1)   and  using  Eqs.    (B-7)    and   (B-14)    in   the  result   it   is   readily  verified 
that  the   coefficients  b    .    for  a  P,   approximation  are 

2    1^2 
bgi   =    1,    bj^   '   -    mi   *    (1-c)    *    (1-c)^}  .   b^j   =  0    . 


(IS).  Using  1-qs.  (B-6)  and  (lJ-3)  in  tq.  (B-1),  defining  R^j  =  A,  U^^  =  B, 
D  for  simplicity  and  setting  to  zero  the  coefficients  of 
and  i(>,i5iti,  (where  the  arguments  have  been  omitted  to  shorten 

the  expressions)  yields  four  nonlinear  equations  for  the  four  unknomis,  A 

through  L): 

1  5    ,.        25      2        147     2   _  ,„    ,., 

T  -  31 "  -  7  '^  *  i  ^^  4  "'^  ^  3i  '^  -  T  '^^  -  m  ^«  =  °  •         e^-i"^ 
-  J  -  ^^  *  7 ' '  ^'^^  -  h^^  *  ^ ''  *  7 ^^  -  m ^"^  - "  •       ^'-''^ 

By  appropriately  manipulating  Eqs.  (B-15)  tlirough  (B-18)  it  can  be  shown  that 
they  are  equivalent  to 

16  -  40C  *  lOOC"  -  147A^  =  0  ,  (B-l'J) 

20CU  -  70AC  +  4aA  -  16D  =  0  ,  (B-20) 

1  -  2C  +  3BC  -  3AU  =  0  ,  (B-21) 

2       2 
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At  this  point  the  P.  variational  boundary  conditions  will  be  considered. 
As  with  tlie  P^  case  one  extra  equation  will  be  required.  This  extra  equation 
again  comes  directly  from  the  spherical  harmonics  equations.  It  is  found  to 
be 

♦  4(x)  +  §y4'2^x)  *   ^  (l-c)4,g(x)  =  0  .  {B-23) 

Using  Eqs.  {B-6),  (U-3)  and  (B-23)  in  liq.  (B-1),  again  defining  Rgj  =  A, 
R.j  =  B,  U.,  =  C,  and  R.,  =  U  for  simplicity  and  setting  to  zero  the  coeffi- 
cients of  't',,S'(i,,,  't'n**]'  *i'^*ii»  ""''  'l'i'5*i»  yields,  after  considerable  algebraic 
manipulation,  four  nonlinear  equations  for  the  four  unknowns,  A  through  D: 
1    I    ,,      .        1  ,,   ,2    r    35   245  ,,   ,  1  „  ^  1225  ^.2 

1    r    35   245  ,,   ,  1  ,,   1  „   1225  ._  ^  161  „.  ^  7  , 

-  -^  AD  -  jil  AU  *  ^  (1-c)  B  =  0  ,  {B-25) 

1    I    35   245  ,,   ,  1  „   1  „   1225  „,,   161  _„  ^  7  , 

*  j^  AD  -  ^  AB  -  j^.  (1-c)  B  =  0  ,  (B-26) 

9    7  „   1225  ,/Z        147  ,,2   „  ,„  ,_, 

-  B-  *  T^  "  *  TtT  "  -  128  '^  -  °  ■  (''-"^ 

By  appropriately  manipulating  Eqs.  (B-24)  through  (B-27)  they  can  be  shown  to 

be  equivalent  to 

144  +  32(l-c)  +  144(l-c)^  ♦  (-280  ♦  980(l-c) )C 

*   2450C^  -  1323A"  =  0  ,  (B-28} 

{22540  -  490(1-c)}c:D  -  3O870AC  +  {35721  ♦  8064(l-c)}A 

-  (23184  +  i)632(l-c)  +  75509(l-c)^)D  =  0  .  (B-29) 
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27  -  54C  *   161BC  -  161AL)  +  32Cl-c)B  =  0  ,  (B-30) 

-  648  +  504li  +  2450U"  -  1323B^  =  0  .  (B-31) 

Equations  (B-28)  tlirougli  (B-31)  are  in  a  convenient  form  for  numerical 
solution.  Given  a  particular  value  of  c,  if  a  value  is  assumed  for  C,  Eq. 
(B-28)  can  be  solved  for  A,  then  Eq.  (B-29)  can  be  solved  for  U,  and  finally 
Eq.  (B-30)  can  be  solved  for  B.  Equation  (B-31)  is  then  used  as  a  consistency 
check.  Equations  (B-19)  through  (B-23)  can  be  solved  in  exactly  the  same  man- 
ner. A  computer  program  is  used  to  compute  the  values  of  A,  B,  C,  and  D.   In 
iterating  to  the  value  of  C  which  best  satisfies  the  consistency  check  equa- 
tion a  relative  accuracy  of  lO"   is  employed.  Since  only  eight  or  nine  digits 
are  retained  in  the  final  results,  this  accuracy  is  deemed  satisfactory.  A 
list  of  the  calculated  values  of  the  constants  A,  B,  C,  and  D  appear  in  Table 
B-I.  The  CMEAM  column  indicates  the  value  of  c  used  to  find  the  corresponding 
constants  A,  B,  C,  and  D. 

The  P.  Marshak  boundary  conditions  when  put  into  the  format  of  Eq.  (B-6) 
yield  the  following  values  for  A,  B,  C,   and  D: 
9 


A 


T 


a  =  -  J  (1  +  J  (!-<=)>  . 
C  =  .j^  (36  +  4(l-c))  , 

Tlie  variational  boundary  condition  constants  for  the  P.  approximation  are 
seen  to  bo  sliglitly  larger  in  absolute  value  than  the  comparable  Marshak  bound- 
ary conditions  when  c  =  1.0.  Using  the  P  variational  boundary  condition 
coefficients  A,  B,  C,  and  1)  corresponding  to  c  =  1.0,  the  Milne  problem  extrap- 
olation distance  is  found  to  be  .708554  mean  free  patlis. 


Table   li-I 

CUNSTANTS    FOR    V/VIATIONAL    BQUNDAKY    CONOITICNS 
(l<IG>JT    HAND    ROUNOARY) 
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P  1  APPROXIMATION 
ABC 
.53259050    -.R6892612     .7A494897 


-1.52200413 


P  A  APPROXIMATION 


CMEAN 


.00 

2 

41228043 

-3. 

34191842 

59843785 

-2.64252082 

.05 

2 

35642719 

-3. 

31748672 

56770363 

-2.62488705 

.10 

2 

30118333 

-3. 

29258481 

53741249 

-2.60691839 

.15 

2 

24657098 

-3. 

26720520 

50758036 

-2.58860980 

.20 

2 

19261313 

-3 

24134083 

47822382 

-2.56995656 

.25 

2 

13933361 

-3. 

21498508 

44936005 

-2.55095430 

.30 

2 

.08675713 

-3, 

18813189 

42100682 

-2.53159909 

.35 

2 

.03490930 

-3 

16077590 

39318254 

-2.51188751 

.40 

.98381658 

-3 

13291247 

36590625 

-2.49181671 

.45 

.93350634 

-3 

10453787 

33919758 

-2.47138453 

.50 

.88400678 

-3. 

07564938 

31307676 

-2.45058955 

.55 

.83534696 

-3 

04624546 

28756463 

-2.42943124 

.60 

.78755677 

-3 

01632585 

.26268256 

-2.40791005 

.65 

.74066682 

-2. 

98589179 

23845245 

-2.38602751 

.70 

.69470847 

-2 

95494615 

.21489668 

-2.36378638 

.75 

.64971370 

-2 

92349363 

.19203804 

-2.34119078 

.80 

.60571501 

-2 

H9154096 

.16989967 

-2.31824632 

.81 

,59703770 

-2 

88509116 

.16556049 

-2.31361617 

.82 

.58840181 

-2 

87862178 

.16125124 

-2.30897242 

.83 

.57980762 

-2 

87213292 

.15697212 

-2.30431513 

.84 

.57125538 

-2 

86562467 

.15272331 

-2.29964438 

.85 

.56274537 

-2 

85909711 

.14850500 

-2.29496024 

.86 

.55427785 

-2 

85255035 

.14431738 

-2.29026279 

.87 

.54585309 

-2 

84598450 

.14016064 

-2.28555211 

.88 

.53747135 

-2 

.83939965 

.13603497 

-2.28082828 

.89 

.52913290 

-2 

.83279592 

.13194056 

-2.27609140 

.90 

.52083802 

-2 

82617344 

.12787759 

-2.27134155 

.91 

.51258696 

-2 

.81953231 

.12384627 

-2.26657883 

.92 

.50437999 

-2 

81287266 

.11984676 

-2.26180333 

.93 

.49621739 

-2 

80619463 

.1158792  7 

-2.25701515 

.94 

.48809942 

-2 

.79949834 

.11194398 

-2.25221440 

.95 

.48002634 

-2 

.79278394 

.10804108 

-2.24740119 

CMEAN 
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.96 

,47199843 

-2, 

,78605158 

.10417076 

-2. 

,24257562 

.97 

,46401594 

-2, 

.77930139 

.10033320 

-2. 

,23773782 

.98 

,45607914 

-2. 

.77253353 

.09652860 

-2. 

,23288789 

.99 

,44818831 

-2. 

.76574817 

.09275714 

-2. 

,22802596 

.00 

,44034369 

-2.75894546 

.08901900 

-2. 

,22315215 

.01 

,43254557 

-2, 

.75212557 

.08531437 

-2, 

,21826660 

.02 

,42479419 

-2. 

.74528868 

.08164344 

-2. 

,21336943 

.03 

,41708983 

-2. 

.73843496 

.07800638 

-2. 

,20846079 

.O'l 

,40943274 

-2. 

.73156461 

.07440339 

-2. 

,20354080 

.05 

,40182318 

-2. 

.72467780 

.07083465 

-2. 

,19860962 

.06 

,39426141 

-2. 

.71777474 

.06730033 

-2. 

,19366739 

.07 

,38674770 

-2. 

.71085562 

.06380062 

-2. 

,18871427 

.08 

,37928229 

-2, 

.70392065 

.06033570 

-2. 

,18375040 

.09 

,37186545 

-2. 

.69697005 

.05690574 

-2. 

,  17877596 

.10 

,36449742 

-2. 

.69000401 

.05351093 

-2. 

,  17379109 

.11 

,35717846 

-2. 

.60302278 

.05015143 

-2. 

,16879597 

.12 

,34990882 

-2, 

.67602657 

,0468274  3 

-2. 

,16379077 

.13 

,342688  75 

-2. 

,66901562 

.04353910 

-2. 

,15877567 

.14 

,33551849 

-2. 

,66199016 

.04028660 

-2. 

,15375084 

.15 

,32839830 

-2. 

.65495045 

.03  707012 

-2. 

,14871647 

.16 

,32132841 

-2. 

.64789673 

.03388981 

-2. 

,14367274 

.17 

,31430906 

-2. 

.64082925 

.03074585 

-2. 

,13861984 

.18 

,30734050 

-2. 

.63374829 

.02763840 

-2. 

, 13355798 

.19 

,30042296 

-2. 

,62665410 

.02456763 

-2, 

.12848735 

.20 

,29355667 

-2. 

,61954695 

.02153370 

-2. 

.12340815 

.25 

,26000204 

-2. 

,58382691 

.00692210 

-2, 

.09789111 

.30 

,22776225 

-2. 

,54782731 

.99325449 

-2. 

.07219285 

.35 

,19686312 

-2. 

,51158873 

.98054847 

-2. 

.04634346 

.'.O 

,16732806 

-2. 

,47515527 

.96881987 

-2, 

.02037550 

.45 

,13917774 

-2. 

,43857439 

.95808245 

-1 , 

.99432388 

.50 

,11242969 

-2. 

,40189670 

.94834763 

-1 . 

.96822571 

.55 

,08709796 

-2. 

,36517566 

.93962430 

-1 , 

.94212008 

.60 

,06319282 

-2. 

,32846719 

.93191850 

-1, 

.91604772 

.65 

,04072041 

-2, 

,29182922 

.92523332 

-1, 

,89005070 

.70 

,01968254 

-2, 

.25532115 

.91956862 

-I, 

,86417200 

.75 

,00007647 

-2. 

.21900323 

.91492102 

-1 , 

,83845509 

.80 

.98189480 

-2, 

, 18293593 

.91128374 

-I. 

.81294343 

.85 

.96512544 

-2. 

. 14717924 

.90864659 

-I , 

.78767998 

.90 

,94975159 

-2. 

.11179197 

.90699607 

-1. 

,76270668 

.95 

.93575191 

-2. 

.07683103 

.90631539 

-I . 

,73806394 

2. 

.00 

.92310068 

-2. 

.04235081 

.90658468 

-I. 

,71379019 
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limployini;  '-I-  C-S)  in  Kq.  (B-O)  it  can  be  shown  that  for  the  P  and  P. 
approximations  tlie  coefficients  b  .  for  tlie  variational  boundary  conditions 


b^,^  =  A.  bj,  =  m.   b^j  =  0.  b3j  =  1/7,  b^j  =  0 
%2  =  '^^  "12 
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APPENDIX  C 


CO  Description  and  Explanation  of  Computer  Programs 


The  group  of  programs  discussed  here  calculate  the  critical  radius  and 
subsidiary  quantities  using  a  P.  approximation  with  any  one  of  the  appropriate 
boundary  conditions  applied  at  the  vacuum  interface.   The  programs  are  written 
for  the  IBM  1410  computer  in  the  FORTRAN  II  language.   Tlie  program  is  set  up 
to  handle  P,  through  P.„  approximations  for  a  bare  spherical  reactor  and  P^ 
through  P  approximations  for  a  reflected  spherical  reactor.   Since  the  size 
of  the  core  storage  of  the  available  computer  is  only  40,000  characters  it  is 
necessary  to  divide  the  complex  program  into  six  phases.   In  addition  each 
phase  is  subdivided  into  a  number  of  subprograms.  The  following  list  indi- 
cates the  order  in  whicli  the  subprograms  are  arranged  within  the  particular 
phases : 


BOLTZMANNl 

110LTZMANN2 

B0LTZMANN3 

INPUT 

MARSMK 

CRITEQ 

POLYCO 

BCMAI« 

SETUPA 

POLYNO 

MARKBC 

CRAM 

EIGEN 

PNPl 

DET 

SETUPG 

VARIBC 

C 

ROOT 

FOURUQ 

CSER 

P  (short) 

NONLIN 

Q 

FACT 

P  (long) 

ROOT 

FACT 

P  (short) 
FACT 

B0LTZMANN4  • 

B0LTZMANN5 

U0LTZMANN6 

CRAM 

INTOUT 

OUTPUT 

SOLVE 

PHIL 

PLOT 

RESIDU 

C 
CSER 

Q 

P  (long) 

FACT 

Each  of  the  programs  in  this  list  will  be  considered  later  in  this  appendix. 
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The  input  data  is  divided  into  two  parts.   The  proper  sequence  for  load- 
ing the  first  part  of  the  input  data  is  shown  in  the  BOLTZMANNl  program  and 
the  INPUT  subproqram.   The  meaning  of  the  symbols  used  in  these  programs  is 
given  in  Table  C-I.  The  second  part  of  the  necessary  input  data  is  discussed 
under  tlie  PLOT  subprogram. 


Table  C-I 


Symbol 


Input  Data 


l3xplanation 


NCASES 
NOIIUUK 
NRBG 


NBC 


ACCU 
IllOLI) 


C(l) 
SIGMA(l) 

C(2) 
SIGMA(2) 

KEFLTU 


Number  of  cases  to  be  considered  in  this  computer  run 

Order  of  approximation,  L 

Number  of  regions  to  be  considered 

=  1  for  bare  reactors 

=  2  for  reflected  reactors  including  the  infinite  black 

reflector  case 

Code  number  for  vacuum-interface  boundary  condition 

to  be  applied 

=   1   for  Marsiiak's  boundary  conditions 

=   2   for  Mark's  boundary  conditions 

=   3   for  variational  boundary  conditions 

=  4   for  infinite  black  reflector  boundary  conditions 

Relative   accuracy  to  be  employed  in  iterations 

=  0  skip  extra  results  and   graohs  used   in  program  testing 

=   1  print  extra  results  used  in  program  testing 

=  0  skip  graphs  used   in  program  testing 

=   1  print   graphs  used  in  program  testing 

Mean  number  of  secondaries,   c,   in  tlie  core  region 

Total  macroscopic  cross   section,   I,   in  the  core  region 

in  cm  _ 

Mean  number  of  secondaries,  c,  in  the  reflector  region 

Total  macroscopic  cross  section,  T,   in  the  reflector 

region  in  cm 

Reflector  thickness,  R  (cm) 


A  sample  page  of  output  is  shown  in  Table  C-II.   Most  of  the  output  is 
self  explanatory.   The  output  shown  in  Table  C-II  is  punclied  on  cards  as  well 
as  being  printed.   Since  output  occurs  in  the  first  and  last  phases  a  special 
card  consisting  of  tiie  number  of  the  case  in  columns  1-5  and  periods  in  the 
remaining  75  columns  is  punclied  for  convenient  separation  of  tlie  cards 
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Table  C-II.  Sample  page  of  output 

P  ">  APPRQXlMAriON 
OF  THE  ONE  VELCICITY  BOLTZMANN  TRANSPOKT  EQUATtON  BY 
THE  SPHERICAL  HAHMCNICS  METHOD  IN  SPHERICAL  r.EOMETRY. 

BAKE  CORE 

VARIATIONAL  V«CUUM  INTERFACE  BOUNDARY  CCNDITIONS 
IN  THE  CURE  C  =  1.0500  ,  SIGMA  =   .1000  /CM. 
THE  ACCURACY  USED  IS  l.OE-07 

THE  COMPUTED  CRITICAL  RADIUS  IS   72.76'.'«  CM. 

CR,    7.276^1  MEAN  FREE  PATHS. 
THE  INTEGRATEP  INWARD  ANGULAR  FLUX  AT  THE 

VACUUM  INTERFACE  IS    7.010E-04 


IN  THE  CORE 


EIGEN  VALUES 
2.53185E  00 


5.171 nE-oi 


GtLiK)  MATRIX 
l.OOOOOE  00 
l.OOOCOE  00 


l.265'»2E-0l   -l.y22B7E-02 
2.58566E-02   -5.20057E-01 


-3.25'.A9E-03 
A.3099'.E-01 


BOUNDARY  CONDITION  MATRIX 

5.32590E-01   -2 .  896'VOE-Ol    O.OOOOOE-00    1.42857E-01 
7.4'i948E-01   -5.0733JE-O1    2.00000E-01    O.OOOOOE-00 

NORMALIZED  CQEFFICTCNTS 

7.95775E-02   -1  .gS'idyE-QS 


DATA  FOR  ANGULAR  FLU>  PLOT  AT  THE  OUTER  BOUNDARY 


MU 

PHKMU) 

MU 

PHI (MU) 

MU 

PHKMU) 

1.0 

2.007E-02 

.3 

8.797E-03 

-.'t 

2.544E-0'i 

.9 

l.85aE-02 

.2 

7.227E-03 

-.5 

-2.885E-0'i 

.8 

1.702E-02 

.1 

5.736E-03 

-.6 

-5.940E-04 

.7 

1.5'.0E-02 

.0 

'..347E-03 

.   -.7 

-6.393E-04 

.6 

1.37'.E-02 

-.1 

3.081E-03 

-.8 

-4.016E-0't 

.5 

1.208E-02 

-.2 

1.962E-03 

-.9 

l.'tl7E-0'. 

.4 

1.0'i2C-02 

-.3 

1.012E-03 

-1.0 

1.013E-03 
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belongini;  to  different   cases.      In  Table  C-II   the   [-IGEN  VALUES  are  the   roots   \^ 
arranjjed  so   that   readinR   from  riglit  to   left  on  each  successive   line   the   first 
value  encountered  is  tlio    imaginary  root    (for  a  core  region)    and  then  the 
remaining  real   roots   in  order  of  decreasing  magnitude.     The  G(L,K),   ^j,(^i,)> 
and  BOUiNUARY  CONUITION,   b    .,   matrices   are  arranged  so  that  the  second  subscript 
is  the   row  number  and  i.  is  the  column  number.      If  H  exceeds  4  then  tlie  remain- 
der of  a  particular  row  appears  on  succeeding   lines.     The  NORI-IALIZED  COEFFI- 
CIENTS are  the  A.  's   and  correspond  in  position  to  the  roots   X,    listed  under 
EIGEN  VALUES.     The   PIIICHU)    column   in  the  inward   flux  section   is  whichever  of 
/(K,u)   or  /(R+U   ,u)    is  representative  of  the  vacuum  interface. 

Approximately  ten  minutes   is   required  to  compute  the   results   shown   in 
Table   C-II   if  this  problem  is   run  with  5   or  more   cases. 

The   sense  switches  do  not   alter  the  program  when  they  are   in  the  off  posi- 
tion.     The  clianges  which  occur  when   they  are   in  the  on  position  are   shown  in 
Table  C-III. 


Switch 


I'hase 


Table   C-III 
Sense  Switches 

Ope  rat  ion  when  Switch  is  on 


liOLTZMANNl         Prints   convergence   in  an   iteration   loop  for 

testing  purposes 
liOLTZMAi>tN2         I'rints   convergence   in  an   iteration   loop   for 

testing  purposes 
liOLTZMANN3         Prints   convergence   in   an  iteration   loop   for 

testing  purposes 
bOLTZMANNS         Sets   lUOLU  in  Table  C-I   to  the  value  of  I 
aOLTZMA'JNb  Sots    III0LD2    in  Table   C-I   to   the   value  of  1 


C.l    UOLTZMANiMl    Program 
This  program  is   tho   control   program   for  the   first   phase.      In  the   first 
phase   tlie   input   data   is    read  and    tlie  preliminary  output    is   printed.      The   roots 
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X.  ,  the  ti,  (A.)  matrix  and  an  initial  lower  bound  estimate  of  the  critical 
radius  are  computed.  The  lower  bound  estimate  of  the  critical  radius  is 


C.2  INl'lJT  Subprogram 
This  subprogram  reads  all  of  the  input  for  a  particular  case  to  be  run 
and  prints  the  preliminary  output  data. 

CT  POLYCO  Subprogram 
Given  a  value  of  c  this  subprogram  sots  up  the  coefficients  of  the  alter- 
nate powers  of  X^   in  liq.  (IS).  Using  Eq.  (14),  Eq.  (15)  can  be  written  as 

\*i^\^  -  (-'''^'ft.i'^'  -<=vw>  =  °  ('-" 

The   Legondre  polynomials   P   (x),   and  the  non-singular  part  of  the   Legendre  poly- 
nomials of  the   second  kind,   W     ,(x),    can  bo   written  as  simple   summations    (20) 

,.    (,)    =    '"f ^         t-^)'t-"-^^J'        x"-2^  (C-2) 

"  i=0     2"i!(n-i)l(n-2i)! 

'V-lf''^    =   J^I^-l'^^'n-.^"'    •  ^'-'^ 

Using  iiqs.    (C-2)    and    (C-3)    in  liq.    (C-1)    it    i;  found  that  the  coefficient  of 
x!'"^'"  wliere  n  is   L+1   and  m  is  an   index  running  from  0  to   [1./2]    in  value  is 

(-1)"    I    (-l)"'(2n-2m)l 
-n     'm!(n-m)! (n-2m) I 

•,      ?        V        V  (-l)J*\2lt-23-2)l(28.-2k-2)l  , 

-  ^'=,i,    .i„  J„it(il,-l-.i)!(l-l-2J)!klj!(n-t-k)!(n-«.-2k)T  "m.j+k' 

where      M  =  the  lesser  of  (li-l)/2  and  m, 
N  =  the  lesser  of  (k-ll)/2  and  m. 
This  is  the  expression  used  ui  this  subprogram  to  compute  the  coefficient  of 


\.        .     Tlie  coefficients   are  arranRed   for  ilescendiriK  alternate  powers  of  \. 
in  the  colunm  matrix  POLY. 


C.4   FOLYNO  Subprogram 
Given  a  value  of  X.    tliis   subprogram  computes  the  value  of  '"i  +  iC^i,)   after 
POLYCO  has  set  up  tlie  coefficients  of  X     in   I'OLY.     This   subprogram  is  used  as 
an  auxilliary  subprogram  by  liKJliN   in   finding  the  roots  of  Eq,    (C-1). 

C.S   lilCJKN   Subprogram 
This   subprogram  supervises  the  finding  of  tlie   roots   \,    of  Eq.    (C-1).      For 
a  core  region  the  signs  of  alternate  coefficients  in  POLY  are  changed  to  find 
the  imaginary  root.     After  the   imaginary  root   has  been  found  the  signs  of  the 
coefficients  in  POLY  are  returned  to  their  original  state  and  the   remaining 
real  roots  are   found.      I'or  a  reflector  region  the   subprogram  finds   the   real 
roots   immediately.      The   real   roots  are  stored  in  order  of  decreasing  magnitude. 
All   roots   are  stored  in  EIC.HiNV. 

C,6  SBTUPG  Subprogram 
This  subprogram  sets  up  the  G   (X  )    matrix.     The  definition  that   G   (A.)   =   1 
for  all   A.    is  used   for  tlie   first   row  of  the  matrix.      By  setting  t  =   1   in  Eq. 
(14)   and  noting  that  W   (A  )   =   1   it   is   readily  apparent  that 
Gj(A^)    =    (c-l)A^    . 

This   last  equation  is  used  to   find  the   G   (A   )'s.     The   remaining  Gj^(Aj^)'s   are 
found  by  applying  the   recursion   relationship   for  these   functions,   Eq.    (13). 

C.7   KUUT  Subprogram 
This  subprogram  is  used  to   find  the  zero's  of  a  function  ^(x).     The 
function   ^(x)    is   specified  by  a  subprogram  whose  name   is   given  as   the   first 
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arBument  of  the  calling  statement.  Given  two  limits  between  whicli  a  zero  of 
^(x)  exists,  tl\is  subprogram  finds  that  zero  to  witliin  the  relative  accuracy 
specified  in  the  calling  statement.  The  method  used  in  finding  the  zero's  of 
^(x)  is  to  apply  a  first  order  iJewton  interpolation  formula  (22).  A  complete 
description  of  this  subprogram  is  on  file  at  the  K.S.U.  Computing  Center. 

C.8  1>  (short)  Subprogram 

This  subprogram  computes  the  value  of  P^CO)-  Setting  x  =  0  in  Eq.  (C-2) 

the  equation  for  F  (0)  is  found  to  be 

0  ,  n  odd 
-„(")  =  <(-l)"/^n!  ^'-'^ 

2"((n/2)!}-  '  ""'"'"   ■ 

C.9   l-'ACT  Subprogram 
Given  a  positive   integer  argument,  n,   this   subprogram  finds  the  value  of 
"n!".      A  detailed  description  of  this   subprogram  is  on   file   in  the  K.S.U. 
Computing  Center. 

C.IU   liULTZMANiN2   Program 

Tills  program  is  the   control   program  for  the   second  phase.     The  second 

phase  sets  up  the   vacuum  interface   boundary  condition  matrix  whose  elements 

are   b    .    in   the   computer   variable   BC. 
i) 

C.ll   MAKSllK  Subprogram 
This   subprogram  in   conjunction  witli  BCMARS  sets  up  the  coefficients  b^^. 
in   liC  corresponding  to  Marshak's  vacuum-interface  boundary  conditions. 

C.12    BCMAKS   Subprogram 
Tliis   subprogram  computes  the  numerical   values   of  b^^.    for  Marshak's  vacuum- 
interface   boundary   conditions   for  wliich 
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(p)U^^'^dU    .  (Cr5) 


EraployiiiR  liq.    (C-2J    in  liq.    (C-5)    and  doing  the   indicated  integration   analyti- 
cally Eq.    (C-5)    becomes 

b    .    =    (-D^^l^'f ' L-iiW2iO^^  .  fC-6) 

"•^  k=0  2\!(i-k:i!(il-2k)!(H+2j-2k) 

Equation  (C-b)  is  used  to  evaluate  the  b  .'s  for  Marshak's  boundary  conditions 

except  when  E,  =  0  in  whicii  case 

b  .  =  l/2j 
oj 

is  used. 


CIS  MAKKBC  Subprogram 
This  subprogram  in  conjunction  with  I'MI'l  sets  up  the  coefficients  b   in 
BC  corresponding  to  Mark's  vacuum-interface  boundary  conditions.   First  the 
roots  of 

are  found  and  then  the  value  of  b  .  is  computed  from 

b,j  =  P,Cu.)  .  (C-8) 

Equation  (C-8)  is  used  to  evaluate  the  b   coefficients,  which  are  stored  in 

*  J 

BC,   for  Mark's  boundary  conditions. 

CIS   VMUBC  Subprogram 

Tliis   subprogram   supervises   tliu   assemblage  of  tlie  b    .    coefficients    for  the 

variational    vacuum- interface   boundary  conditions.      Tlie   variational   coefficients 

b    .    are  defined  in  Appendix   u.      I'or  the   P,   and   I',   aoproximations   it   is  neces- 
l)  '  3  4"^ 

sary  to   solve   a   set   of   four  nonlinear  ecjuations.      This   subprogram  sets   up   a 
matrix,   CO,   of  tiie   coefficients   in  these  nonlinear  equations   and  supervises  the 
solving  of  tl\ese  equations. 
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C.16  FOUREQ  Subprogram 
This  subprogram  is  an  auxilliary  subprogram  used  by  VAIUBC  in  finding 
the  solution  of  the  four  nonlinear  equations  for  the  variational  boundary 
conditions. 

C.17  NONLIN  Subprogram 
Given  a  guessed  value  of  C  tliis  subprogram  solves  three  of  the  nonlinear 
equations  representing  the  variational  boundary  conditions  and  uses  the  fourth 
equation  as  a  consistency  check, 

C.18  i'  (long)  Subprogram 

This  subprogram  computes  the  values  of  the  Logendre  polynomial  P„(x). 

In  order  to  obtain  good  accuracy  for  small  values  of  x, Eq.  (C-2)  is  rearranged 

into  the  form  ' 

,  ,,n  m  ,  ,,k,.-  ^  11  ^  I  n-2m-2k 
P  rvi  -  C-i)   V  (-1)  (2n-2m-2k)lx  .og, 

'^n^'''  ■  ■^r~^ijj(m-k)  I  (n-ii.*k) !  (n-2m*2k) !  ^''  ' 

where 

m  =  [n/2]  . 
Equation  (C-'J)  is  used  to  compute  P  (x)  except  when  x  =  0  in  which  case  Eq. 
(C-4)  is  used.  A  detailed  description  of  tliis  subprogram  is  on  file  at  the 
K.S.U.  CompuLing  Center. 

C.l!)  I10LTZMANN3  Program  ' 
This  program  is  the  control  program  for  the  third  phase.  The  third  phase 
sets  up  tlie  T  matrix  in  Eq.  (44)  in  the  computer  variable  A  and  iterates  to 
find  a  critical  radius  which  satisfies  Eq.  (60). 
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C.20  CRiTEQ  Subprogram 
This   subproRram  is  an  auxiliary  subprogram  used  by  B0LTZMANN3   in  iterating 
on  Eq.    (60)    to   find  the   critical   radius. 

C.21    SHTUI'A  Subpro;;ram 
Tins   subprogram  sets  up  the  T  matrix  in  the  two-dimensional  computer 
variable  A  for  the  particular  reactor  system  being  considered. 

C.22   CRAM  Subprogram 
'nils  subprogram  performs  the  first   step  in  a  Crout   reduction  method  of 
solving  a  matrix  equation.      Tlie   first   step  in  the  Crout   reduction  is  to  reduce 
the  given  Matrix  to  an  upper  right   triangular  matrix.     The  Crout   reduction 
formulas   for  this   first   step  arc   given   as   (7J 

a!  .    =  a.  .    -     y  al,  a,'  .  for  i   >  j    , 

where  tlie  primes  denote  the  transformed  elements.  This  subprogram  has  special 
provisions  which  allow  imaginary  row  interchanges  so  as  to  maximize  the  diag- 
onal in  the  reduction  process.   A  detailed  description  of  this  subprogram  is 
on  file  at  the  K.S.U.  Computing  Center. 

C.2.'5  UKT  Subprogram 
'ITiis  subprogram  calculates  tile  determinant  of  a  given  matrix  after  CRAM 
has  reduced  the  matrix  to  an  upper  right  triangular  matrix.  Since  the  matrix 

is  in  the  upper  right  triangular  form  the  determinant  of  a  transformed  n  by  n 

matrix  A'  is 

|A'|  ^  C-D'  n  a! 
i  =  l 
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where  H  is  the  number  of  imaginary  row  interchanges  used.  A  detailed  descrip- 
tion of  this  subprogram  is  on  file  in  the  K.S.U.  Computing  Center. 

C.24  C  Subprogram 
This  subprogram  calculates  the  values  of  the  C  (x)  functions.   For  real 
arguments  Hq.  (A-G)  is  used  whereas  for  imaginary  arguments  whichever  of 
Eqs.  (A-11)  or  (A-12)  is  appropriate  is  used. 

C.25  CSER  Subprogram 
This  subprogram  evaluates  the  finite  sums  in  Eqs.  (A-11]  and  (A-12)  for 
the  C  subprogram, 

C.26  Q  Subnrogram 
Tliis  subprogram  computes  the  values  of  Qj(x)  by  using  Eq.  (A-3).  A 
detailed  description  of  this  subprogram  is  on  file  in  the  K.S.U.  Computing 
Center. 

C.27  U0LTZMA1^IN4  Program 
This  program  is  the  control  program  for  the  fourth  phase.  The  fourth 
pliase  having  been  given  the  critical  radius  by  the  third  pliase  assumes  A  is 
one  and  solves  for  the  remaining  A,  and  A,  constants  in  the  column  matrix  X. 

C.28  SOLVE  Subprogram 
This  subprogram  solves  a  matrix  equation  of  the  form 

AX  =  B  (C-10) 

for  tlie  column  matrix  X  after  CHAM  has  reduced  the  matrix  A  to  an  upper  right 
triangular  matrix.  The  formulas  used  in  this  second  step  of  the  Crout  reduc- 
tion are  (7) 

^1  =   IT-  f"i  -  X'^WO 
11       k=l 
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and 


X.  =  b!  -   y   a!,  X, 
1    1   ,  '■.  ,  ik  k 

where  tlie  primes  denote  transformed  quantities.  A  detailed  description  of 
this  subproRram  is  on  file  in  the  K.S.U.  Computing  Center. 

C.2y  RHSIDU  Subprogram 
This  subprogram  minimizes  the  error  incurred  in  using  the  Crout  reduction 
method  of  solving  a  matrix  equation.  After  CliAM  and  SOLVE  have  been  used  to 
find  a  solution  X  this  subprogram  computes  the  product  of  the  matrix  A  and 
the  solution  X  and  subtracts  the  result  from  the  column  matrix  B  in  Eq.  (C-IO). 
Then  CRAM  and  SOLVE  are  again  used  to  solve  Hq.  (C-10)  for  a  new  column  matrix 
X.  This  result  is  used  to  reduce  the  columji  matrix  B  in  the  same  manner  as 
before  and  then  it  is  added  to  the  old  column  matrix  X.  Then  a  new  column 
matrix  X  is  found  by  solving  Hq.  (C-10).  The  iteration  is  continued  until  the 
column  matrix  X  determined  after  the  n   iteration  is  negligible  (to  within 
the  accuracy  of  the  computer)  with  respect  to  the  sum  of  the  results  of  the 
first  n-1  iterations.  A  detailed  descripticin  of  this  subprogram  is  on  file 
in  the  K.S.U.  Computing  Center. 

C.30  UOLT2MANN5  Program 
This  program  is  the  control  program  for  the  fifth  phase.  The  fifth  phase 
computes  all  of  tlie  data  which  will  be  printed  in  the  final  output.   The 
amount  of  data  which  is  computed  is  governed  by  the  input  constant  lilOLD  or 
sense  switch  1. 
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C.31  IMTOUT  .Subprogram 
This  subprogram  computes  all  of  tlie  necessary  data  for  the  final  output. 
The  amount  of  computation  done  is  controlled  by  the  input  constant  IHOLD  or 
sense  switcli  1. 

C.J2  I'HIL  Subprogram 

This  subprogram  computes  the  value  of  ^^(r).  If  r  represents  a  point 

which  is  in  the  core  region, 

[(Ul)/2] 
^  (r)  =  (24*1)   I       A^(;^(X^)Cj^(!:r/X|^)  (C-U) 

k=l 

whereas  if  r  represents  a  point  which  is  in  the  reflector  region, 

2[(U1)/21_ 
7^(r)  -  (2U1)     I         \(^i(\)Qi(!^r/X|.)  .  (C-12) 

k=l 

The  A*  terms  are  omitted  from  liqs.    (C-11)    and   (C-12)    since  they  arise  only  at 

an  interfacial  boundary  between   two  media. 

C.33    liOLTZMVJNO   Program 
This  program  is  the  control   program   for  the   sixth  phase.     The  sixth 
phase  prints  and  punches   the  computed  results.      In  addition   graphs  of  parts 
of  the  results  may  be  printed  under  tlie  control   of  the  input   constant   IH0LD2 
or  sense  switclt   1. 

C.34   UUTi'UT  Subprogram 
This  subprogram  supervises  the  printing  and  punching  of  the  computed 
results.     Wlien   graphs   are  asked   for  by  having  IH0LI)2   set  to  one  or  sense 
switch   1   on,   three   graphs   are  printed   for  a  bare   spherical   reactor  whereas 
seven   graphs   are   priiited    for  a    reflected   reactor.      For  a  bare   reactor  plots 
are  made  of  tlie  total   neutron   flux   in  tlie   core  as   a  function  of  tlie   spatial 


in   the  core   as    a   fujiction  of  tlio   spatial   variable.      In   addition,    for  a 
reflected  reactor  plots   are  made   of  tlio  total  neutron   flux  in  the   reflector  as 
a  function  of  the  spatial   variable,   tlio  individual   ^   (r)'s   in  the   reflector  as 
a  function  of  the   spatial   variable,    and  of  the   angular  flux  at   the   core- 
reflector  interface  as  viewed   froin  eacli  side  of  the  interface, 

C.35   PLOT  Subprogram 
This   subprogram  plots   a  jiraph  of  some  of  tlic  computed  results   in  tlie 
printed  output.      In   order  to  print  titles   and   axes   of  graphs, control   cards   are 
read  in.      A  list   of  the   control   cards  used  in  tliis  work  are  shown  in  Table 
C-IV.      A  detailed  description   of  this   subprot;ram  is   on    file   in  the   K.S.U. 
Computing  Center, 

C.ib   Special   Macliine    Language  Subprograms 
There   are   tliree   special   machine    language   subprograms  which   are  used  in 
this   program.      The   machine    lanjjuago   subprogram  CtliiCKI    is   a   routine  whicii  checks 
the   error  indicators   and  prints   appropriate    3rror  messages   if  necessary,      liXIT 
is   a  machine    language    routine  wiiich  terminal.cs   program  control    and  gives   con- 
trol  to  the   comjiuter  monitor.      I'Ue   machine    l.'Lnguage   subprogram   INfJUIK  checks 
the   inquiry   reijuest   key  on   the   console   of  tlie   computer  allowing  alteration  of 
the   sense   switcli   settings.      A  detailed  description  of  eacli  of  these   subprograms 
is   on   file   in   the   K.S.U.    Computing  Center. 


93 

■  L ; 


I    H    I  I  I  I  I    "J  i 

+0+       +       +       +       +a+ 
z    .  «    .        .        .  a    .        •  -    • 

O  1-1  I-  (-  LU  l/l 

—         Z  3  OC         U.   3  Q 

1/1      uj  o  a  a  —  oc 

z      £  a£  u.           Di  irt  Q 

Ui  1/)  —  o  uo  z  H 

rOQQ  300LU  u 

«■-.         1-  Z-i">-OC  UJ 

00_J  OQt-        O  -J 

^  a:  —      i<ii-iKZ  <uj< 

<         QE        l-U        Z>3LaJt-Q:i- 

-.j<ou.ujz_iu.ozuu  uj 
o<t£a:oi3<      oct-iiul 

<0  U-  h-U.U<ll.  f-IH 
a  —  UJ         Z  "-  J         I- 

KIU-LJU^t-1/^UJDU-         U- 

aJ«^-o  —  o  HH<ii_)a£a2:a 
X!^       21—       too:       zu-       o 

a>-zz  ai-^zu-Z 

fc  U.11.       i/iZDi-iuju-i-iyia:—       ►-< 

S         oaz»-<u.i/icT£:oo  —  <i/i>i^ 
S.  OQ      oa      ujOQOasa 

Z  Z  "        <OOZ-'        zoco 
>^  aij>-u-         *■         □LLU-Z3*'Q  — 

^  —  — lUOi/)         LU^UJOO  Z 

_  i-H-z      <Dri-a:      co3or> 

u  ounz      xi-o      z      tor 

ifl  z<ii-auj  <iua-i      oD 

3  oa      -.o      za:Xi-.< 

^  o  u.  u  o       < 

{;  <       i/><a;       z       z<<       1-1 

u  I/)           u.(-""           u.o:< 

^  <a;zi-zuju. 

7?  a      —      scrocri-a 

^  ujh-                cnco^-'oszuj 

«  a       o       £       — r-i-r~'-       >- 

C  OUJ3oi'C^<l  z 

o         o      -J      3      h-i/>inir\uj      « 
'-'  u        c        u    <l  —  1  i. 

LU       LU       <       i-imoi'^^H-       ili 

I  OC          >         Q(M(-(N  I 

>  t-  _  1/1  -J  H-  H- 

►-'  •LU*LtJ*h-0»-'0<*           • 

A  Z  X          I         Z         Q  I- 

"  „  I-         I-         UJ  Z  < 

4)  E        I-  □ 

'-i  ZZI-OZ  —  Z 

-9         o'-'<xuj>-a 

2         -  r       3      - 

■^         t-      z      X      _i(-a>-cD      I- 

oxax3      <zi:z'—      3 

CQ=>~Z>_I         -'UJ         LUo;         CD 

—  _iH_iu.xi-r_ji:t-x  —  x 
au.3u-  3<o<oi/i3oi:3 
I-      £c      a:_io-£  —  s:~_ii-_i 

l/)QH-a<ll.l/l  h-  PU.l/)U. 
^UJQCUJJ  j<_j  -* 

Ois^t—  f^3o£-J<Q-<xci:QQ: 
•-"i/>>-'0<<3i/l3D<        < 
X_|i-i_IZ_J3Q         Q-IJK-J 

_jz:      X      o-'><>      o_jcj> 

u-ocxa       z>  —  oi-ioczu-z 

ODD      <^aQa<<       < 

_]Z_I2  Q2  —  Z_l         a 

<  U.  Z»-'>"3         < 

t_  -*        —        O         -J 

o       -I  D       z      n 

»-      <  z      <      o 


C.37  Computer  Program  Logic  Diagram 


94 


rREAD  PHASE 
FROM  TAPE  5 


SET 
ICASE :  I 


READ  NCASES 
FROM  TAPE 


;es\  I 


READ     COMMON 

AREA 
FROM  TAPE   4 


X 


SET   UP 

APPROPRIATE 

BOUNDARY 

CONDITION   MATRIX 


CONTINUED 

ON 
NEXT   PAGE 


TRANSFER  CONTENTS' 
OF  TAPE  8 
TO  TAPE  4 


SET 
IREG  =  NREi^ 


ZE 


CALL    POLYCO 
CALL    EI6EN 
CALL    SETUPG 


SET 
IREG  =  I 

~T~ 


MOVE    GC    AND 
EIGENC    MATRICES 
TO  G   AND   EIGENV 


TFIANSFER   CONTENTS 
OF  TAPE    8 
TO  TAPE  4 


/READ  PHASE 
\FROM  TAPE  i>J 

C  REWIND 
'ES  4  AND  8 


* 

SET  UP  MATRIX 

A 

AND  FIND 

CRITICAL    RADIUS 

JlL 


READ    COMMON 

AREA 
FROM  TAPE  4 


COMPUTE 

ESTIMATE  OF 

MINIMUM 

CRITICAL     RADIUS 


f         READ 

NCASES 
VFROM  TAPE  4 


SET 

KCASEsKCASEtl 


SET 
KCASE  =   0 


95 


CONTINUED 

FROM 

PREVIOUS    PAGE 

/read  PHASE  4 
\FROM   TAPE  5 


96 


C.38  Computer  Program  Listing 

HOP       BOLTZMANNl  1  <. 

C  THIS    PRUGRAM     IS    THE    CONTROL    PROGRAM    FOR    PHASE    1 

DIMENSION    P0LY16)  ,  E  I  GENV  t  "i  )  ,  n  (  1 1  ,  5)  ,  GC  I  1 1 ,  5  )  ,C  (  2  )  ,  SI  GMA(  2  1  ,E  IGENC  ( 
15) 
DIMENSION    BC( 11. 5) 

COMMON    POLY,N,EIGENV,G,NREG,NQC,C,SIGMA,REFLTH,EIGENC,GC,ACCUR,HOL 
COMMON    H0L2,0C 
REW IND4 
READ    INPUT    TAPE    5,6,NCASES 

6  FORMAT! I^) 

WRITE  TAPE  A.NCASES 
DO  1001  ICASE=1.NCASES 
TYPE  7,ICASE 

7  FORMAT(  16X1'.HPR0DLEM    NUMBER,  15) 
WRITE    OUTPUT    TAPE    ft.lUl.ICASE 

nil    FORMAT!  2H1$,  I5,75H 

1 ) 

C      INITIALIZATION  OF  THE  COMMON  AREA 
1PR0B=1 
NCV2=2 
NPl  =  5 

POLYd  )=0.0 
REFLTH  =  O.C) 
DO  8  K=l,2 
C(K)=0.0 

8  SIGMA(K)=0.0 
DO  9  K=1,N0V2 
E1GENV(K)=0.0 
EIGENC(K)=0.0 
PQLYIKtl)=0.0 
DO  9  L=l,NPl 
G(L,K)=O.U 
GC(L,K)=0.0 
BC(L,K)=0.0 
IREG=N0V2+K 

9  GIL, IREG)=0.0 
CALL  INPUT 

C      THIS  SECTION  OF  THE  PROGRAM  "^ETS  UP  G(L,K)  MATRIX  AND  FINDS  THE 
C        ROUTS  LAMBDA-SUB  K. 

00  19  KKEG=l,NREG 

IREG=NREG*1-KREG 

IF(C( IREG)-1. ) 12, 10, 12 

10  WRITE  OUTPUT  TAPE  6, 11 

11  F0RMAT(5XJ2HC  =  1, PROBLEM  CANNOT  BE  EXECUTED) 
IPRUB=2 

GC  TO  1000 

12  IF(N-2«(N/2)  )16,  13,  16 

IJ  1F(C( IREG)-(FL0AT(N+1)«PIN,0.) )»«2)16, 14,14 

14  WRITE  OUTPUT  TAPE  6  ,  15  ,  C  (  IRE*" )  ,  N 

15  FCRMAr(5X16HrHE  VALUE  OF  C  = , IPE I  5. 7 , 3X9HIN  THE  P  ,I3,14H  APPROXIM 


97 


lATILiN,  IXAbHlS    GREATER    THAMQ'^    EQUAL    TO    (  (  Nt  1)  •?  (  N,0)  )  »»2) 
IPRGB=2 
GO    TO    1000 

16  CALL  PQLYC0(N+1,C( IREG) iPOLY) 
CALLEIGEN(C( IREG) ) 

CALL    SETUPG(C( IREG) ) 
IF(  IREG-n  19,19,  17 

17  NDV2=lNtl)/2 
NP1=N+1 

aO    18  K=1,N0V2 
EIGENV(K)=EIGENC(K) 
00  18  L=1,NP1 

18  G(L,K)=GC(L,K) 

19  CONTINUE 
POLY( l)=HUL 
P0LYI2)=H0L2 

C      COMPUTE  AN  ESTIMATE  OF  THE  MINIMUM  POSSIBLE  CRITICAL  RADIUS  FROM 
C        A  P-2  CALCULATION. 

ALAMBD=SQRT( (9.-^.»C(l) )/(15.»(C( 1)-1.) ) ) 

K-1.E*05 

SQRT3=SQRT(J. ) 

20  RLAST=R 

R  =  ATAN( l./ll./R-l./(    SQRT3«tC(  1  )- 1 . ) »ALAMBD) ) ) 
IF(  1030,  30,21 
30  R  =  Rt3.1^1'3927 

IFI ABSI (R-RLAST)/R)-ACCUR  121,20,20 

21  GO  TO  (22,23),NKEG 

22  R=.9«R 
GO  TO  2^1 

23  R=.3»R 

24  R=R«ALAMBO/SIGMA( I ) 

1000  CALL  CHECK! 

1001  WRITE  TAPE  A , I PROB , POL Y , N .E IHENV , G,NREG, NBC , C , S IGMA , REFLTH.E IGENC, 
1GC,ACCUR,HC,R 

REWIND  4 
CALL  EXIT 
STOP 
END 

BOP   INPUT  14 

SUBROUTINE  INPUT 
C      THIS  SUBPROGRAM  READS  IN  THE  INPUT  DATA 

DIMENSION  POLY (6) ,E1GENV( 5 ) , G I  1 1 , 5 ) , GC t 1 1 , 5 )  ,C ( 2 ) , SI GMA( 2 ) t E IGENC ( 
15) 
COMMON  POLY,NQROER,EIGENV,G,NKEG,NflC,C,SIGMA,KEFLTH,EIGENC,GC,ACCK 
COMMON  H0LD,HQLD2 
READ  INPUT  TAPE!>,2,N0R0ER,NRFG,NBC,ACCR,  IHOLD,  IH0LD2 

2  FCRMAT(3I5,E15. 8,215) 
HOLD=IHULO 
M0LD2=IHOLD2 

READ  INPUT  TAPE  5 , 3, C (  1 ) , S IGMA I  1 ) 

3  FORMAT! 3E15. 8) 
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CO  TU  (S.AI.NREG 
A  1F(NRC-4)'.1,42,41 
'»2  C(2)=0.0 

SlGMAt2)=SIGMAI I ) 
GO  TO  5 
C      THE  CONSTANTS  FOR  THE  REFLECTOR  ARE  NEEDED  ONLY  IF  THIS  IS  A 
C      REFLECTED  REACTOR  CASE. 

41  READ  INPUT  TAPE  5 , 3, C ( 2 ) , SIGWA ( 2 ) .REFLTH 
C      THIS  SECTION  PRINTS  A  STATEMENT  OF  THE  CASE  BEING  CONSIDERED. 

5  WRITE  OUTPUT  TAPE  6,6,NURDER 

6  FORMAT! 3H-$1 , 3 IX IHP , I  2, IX 13HAPPR0XIMAT I  ON) 
WRITE  OUTPUT  TAPE  6,7 

7  FORMAT! 3H  $  ,15X51H0F  THE  ONF  VELOCITY  BOLTZMANN  TRANSPORT  EQUATIU 
IN  BYl 

WRITE  OUTPUT  TAPE  6,71 

71  FORMAT! 3H  $  ,15X53HTHE  SPHERICAL  HARMONICS  METHOD  IN  SPHERICAL  GEO 
IMETRY.I 

IF!NBC-4)72,8,72 

72  GO  TO  18,10),NREG 

8  WRITE  OUTPUT  TAPE  6,9 

9  F0RMAT!3H0$0,20X9H8ARE  CORE  I 
GO  TO  12 

10  WRITE  OUTPUT  TAPE  6,11 

11  F0RMAT!3H0$0,20X14HREFLECTE0  CORE)  , 

12  GO  TO  113, 15, 17, 25), NBC 

13  WRITE  OUTPUT  TAPE  6,14 

14  F0RMAT!3H  $  , 20X44HMARSHAK  VACUUM  INTERFACE  BOUNDARY  CCNDITIGNS) 
GO  TO  19 

15  WRITE  OUTPUT  TAPE  6,  16 

16.F0RMATI3H  $  ,20X42HMARKS  VACUUM  INTERFACE  BOUNDARY  CONniTIONS) 
GC  TO  19 

17  WRITE  OUTPUT  TAPE  6,  18 

18  FURMAT!3H  $  , 20X4BHVAR I  AT  ION«L  VACUUM  INTERFACE  BOUNDARY  CONDITION 
IS) 

GO  TU  19 

25  WRITE  OUTPUT  TAPE  6,26 

26  F0RMATI3H  $  , 20X40H1NF INI TELY  REFLECTED  WITH  A  BLACK  MEDIUM) 

19  WRITE  OUTPUT  TAPE  6, 20, C I  1 ) , S IGMA I  1 ) 

20  FCRMATI3ri  $  ,20X15HIN  THE  COPE  C  =,F7.4,1X9H,  SIGMA  =  ,F7.4 , 1 X4H/CM 
1.) 

GO  TO  !24,21),NREG 

21  IFINBC-4)27,24,27 

27  WRITE  OUTPUT  TAPE  6 , 22 , C ! 2) , S IGMA ! 2 ) 

22  FDRMATI3H  $  ,20X20HIN  THE  REFLECTOR  C  =,F7.4,1X9H,  SIGMA  =,F7.4,1X 
14H/CM. ) 

WRITE  OUTPUT  TAPE  6,23,REFLTH 

23  F0RMAT!3H  $  ,20X26HTHE  REFLECTOR  THICKNESS  1 S , F 7. 2 , 1 X3HCM. ) 

24  WRITE  OUTPUT  TAPE  6,31,ACCR 

31  F0RKATI3H  $  ,20X20HTHE  ACCURACY  USED  IS.lPEB.l) 
RETURN 
END 
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BOP   PQLYCO  I') 

SURKOUTINt  PQLYCO(N,C,PaLY) 
C      THIS  PROGRAM  SETS  UP  THE  COEFFICIENTS  OF  THE  POLYNOMIAL  NECESSARY 
C        FUR  OETFRMINING  THE  VALUES  OF  LAMBDA.  THE  COEFFICIENTS  ARE 
C        STORED  IN  POLYI I ). 

DIMENSION  POLY(l) 

NGV2=N/2+l 

DOB  NM-1,N0V2 

M=NM-l 

POLY(NM)=( (-1. )»*M)»FACT(2»(N-M) )/( FACT ( M) "FACT ( N-M) •F ACT ( N-2»M) 1 

DO  b    L-1,N 

LMlOV2=^lL-l)/2  +  l 

NML0V2=(N-L)/2tl 

!FINM-LM1QV2)1,2,2 

1  LH10V2=NM 

2  00  6  NJ=l,LMlOV2 
J=NJ-l 
IF(NM-NML0V2)3,'.,'. 

3  NMLQV2=NM 

4  DO  6  NK=1,NMLUV2 
K=NK-1 
1F(M-K-JI6,5,6 

5  P0LYINM)  =  P0LYINM)-2.»C»FACT(?»(L-l-J) ) "FACT  I  2* ( N-L-K )  )»(t-l.l»»(K+ 
IJ) )/(  FACT! J)»FACT1L-1-J ) "F AC T I L-1-2* J ) »FAC T t K ) 'FAC T ( N-L-KI 'F ACT  I 
2N-L-2»K)»HL0AT(L ) ) 

6  CONTINUE 

8  P0LY(NMI=POLY(NM)«( |-.5)*«N) 
RETURN 
END 

BOP   POLYNO  lA 

FUNCTION  POLYNQIX) 
C      THIS  SUBPROGRAM  IS  USED  AN  AI.'XILLIARY  SUBPROGRAM  BY  EIGEN  IN 
C      FINDING  THE  ROOTS  LAMBDA-SUB  K. 
C      IT  ACTUALLY  COMPUTES  THE  VALUE  OF  G(L+l,X). 

DIMENSION  P0LY(6) 

COMMON    POLY.N 

NCO=     (N+n/2  +  1 

P0LYN0=POLYlNCO) 

X2=X*X 

FMULT=1. 

UO  I  1=2, NCO 

NSUB=NCO-ltl 

FMULT=FMULT»X2 
1  POLYNO=POLYNO+POLY(NSUBI»FMULT 

RETURN 

END 

HOP   EIGEN  I'. 

SUBROUTINE  EIGEN(C) 
C      THIS  SUBPROGRAM  FINDS  THE  ROOTS  LAMBDA-SUB  K 

DIMENSION  P0LY(6),EIGENV(t)),(;(  ll,'jl,GC(U,5),U(2),SIGMA(2),EIGENC( 
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15) 

COMMON  POlY,N,EIGENC,GC,NREG,NBC,D,SIGMA,REFLTH,EIGENV,G,ACCUR 

NPl=lNtl)/2tl 

IF  THE  VALUE  OF  C  IS  GREATER  THAN  I , F I  NO  THE  IMAGINARY  ROOT  FIRST. 

IF(C-I. )1,2,2 

1  KIMAG=0 

CO  TO  3  ., 

2  KIMAG=1 

00  U  1=1, NPl 

IF( l-2«( 1/2) )ll, 10, II 

10  PGLYl 1 )=-POLY( I) 

11  CONTINUE 

3  J=l 

9  X=1C0. 

ANS2=P0LYN0(X) 

DO  5  1=1,205 

ANS1=ANS2 

IF(  i-90)9'.,9'«,92 
92  IFI  1-180)93,93,91 
91  DELT=.0^ 

GO  TO  95 
9  3  DELT=.l 

GO  TO  95 
9^,  D£LT=1. 
95  X=X-OELT 

ANS2=P0LYN0(X) 

IFISIGNl  1.  ,ANS2)-SIGN(  l.,ANSl  )  )4,5,'i 

POLYNO 
'i    EIGENVlJ)=ROOT(  POLYNO,  X  +  DELT,X,0.0,ACCUR) 

J  =  Jtl 

IF(KiMAG-l )5,6,6 

5  CONTINUE 
81  RETURN 

6  KIMAG=0 

00  8  1=1, NPl 

IF( l-2«( 1/2) )8,7,8 

7  POLYd  )=-POLY(  I) 

8  CONTINUE 

1F( J-NP1)9,81,81 
END 

HOP   SETUPG  l-i 

SUBROUTINE   SETUPG(C) 

THIS  SUBPROGRAM  SETS  UP  THE  ELEMENTS  OF  THE  MATRIX  G(L,K) 
DIMENSION  P0LY(6) , E I GENV I  5 ) , G (  ll,5),GCIll,5),D(2),SIGMA(2)  .EIGENCI 
15) 
COMMON  POLY,N,EIGENC,GC,NREG,NBC,0,SIGMA,REFLTH,EIGENV,G,ACCUR 
N0V2=(N+l)/2 
NP1=N+1 
DO  2  K=1,N0V2 
G(1,K)=1. 
G(2,K)=(C-1. )«EIGENV(K) 
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IF(N-l)2,2,l 

1  DO    20    L=3,NP1 
SIGNA=-1. 
IF(C-1. )2U,18, 18 

la     IFU-1  )21, 21,20 

21     IF( L-l-2«l(L-l)/2) )20, 19,20 

19  SIGNA=l. 

20  G1L,K)  =  (FL0AT(2»L-3)»EIGE1MV(K).G(L-1,K)»SIGNA-FL0AT(L-2)«G(L-2,K)  ) 

1/FL0AT(L-1) 

2  CONTINUE 
i    RETURN 

END 

BOP   ROOT  I'" 

FUNCTION  r<OOT(  DUMMY,  RFB  ,  RFC  ,  RFAE  ,  RFRC  ) 
C      THIS  SUBPROGRAM  COMPUTES  A  ROOT  OF  A  GIVEN  FUNCTION. 
C        DUMMY  IS  THE  NAME  OF  THE  Sl'BPRUGRAM  WHICH  REPRESENTS  THE  FUNC. 
C        RFB  IS  THE  LOWER  LIMIT  OF  THE  RANGE  IN  WHICH  THE  ROOT  IS  TO  BE 
C  FOUND. 

C        RFC  IS  THE  UPPER  LIMIT  OF  THE  SAME  RANGE. 
C        RFAE  IS  THE  ABSOLUTE  ACCURACY  TO  BE  ITERATED  TO. 
C        RFRE  IS  THE  RELATIVE  ACCURACY  TU  BE  ITERATED  TO. 

C        WITH  SENSE  SWITCH  6  ON  THE  CONVERGENCE  AT  EACH  TRIAL  IS  PRINTED. 
JRFS=1 

RFFB=DUMMY(RFB) 
RFFC=DUMMY(RFC) 
I F(RFFC«RFFB 19122, 9122, 9102 
9102  WRITE  OUTPUT  TAPE  6 , 200 , RFB , RFC 
200  FORMAT! 1X76HLIMITS  GIVEN  TO  ROOT  FUNCTION  GENERATE  FUNCTIONAL  VALU 
lES  WITH  THE  SAME  S  IGNS  ,  1 1  H,  L  IM  I  TS  ARE  ,  IPE  18.  T.'iXaHANO  ,F18.  7) 
RUOT=0.0 
RETURN 

9122  RFA=RFC 
RFFA=KFFC 

9123  IF(ABS(RFFB)-ABS(RFFA))913l,913l,9130 

9130  RFC=RFB 
RFR=RFA 
RFA=RFC 
RFFC=RFFB 
RFFB=RFFA 
KFFA=RFFC 

9131  RFD=0.5»IRFB-RFAI 
RFT=RFAE+KFRC«ABSIRFDI 
IF(RFFB)9156, 9135, 9156 

9135  R0OT=RFB 
RETURN 

9156  IF( IKFK)9137, 9157, 9137 

9157  IRFK=3 

9138  RFV-RFD 

GC  TO  91',0  1J 

9137  IF(KFFB-RFFC19139, 9138, 9139  j 

9139  RFV=(RFB-KFCI»RFFB/(RFFB-KFFC) 
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9140  IF(ABSIRFC)-RFT)91'iJ,9l'.3,91'i4 

9143  JRFS=2 

GC  TO  9147 

9144  IF(ABS(HFV)-RFT)9149,9147,9147 

9147  IF( I RFD-RFV)«RFV 19148,9152, 9152 

9148  RFX=RFB-RFO 
GC  TU  9153 

9149  IF(I<FD)9151, 9150, 9150 

9150  RFV=KFT 

GC  ro  9152 

9151  RFV=-RFT 

9152  RFX=RFB-RFV 

9153  GC  TO  (9101,9100 l.JRFS 
9100  RCOI=RFX 

RETURN 
9101   CALL  INQUIR 

IF(SENSE  SWITCH  6)9111,9112 

9111  WRITE  OUTPUT  TAPE  6, 1000 , RFA ,RFB , RFX 

1000  FORMAT! IX45HPRESENT  LIMITS  OF  INTERVAL  IN  ROOT  SUBPROGRAM, IPEIS. 7 , 
14X2HTQ,E18.7,4X16HNEXT  TRY  WILL  BE, CIS. 7) 

9112  RFC=RFB 
RFB=RFX 
RFFC=RFFB 
RFFH=DUMMY(RFX) 
IF(RFFA«RFFB)9158,9131,9159 

9158  IRFK=IRFK-1 
GO  TO  9123 

9159  1RFK=2 

GO  TO  9122 

STOP 

END 

BOP   P  14 

FUNCTION  P(N,X) 
C      THIS  SUBPROGRAM  COMPUTES  THE  VALUE  FOR  AN  LTH  ORDER  LEGENORE 
C        POLYNOMIAL  WITH  A  ZERO  ARGUMENT. 
C      LEGENORE  POLYNOMIAL  FUNCTION  CALCULATOR,  ASSUMING  X=0.n  . 

3  M=N/2 

4  IF(N-2»M)7,8,7 

8  P=(  (-1.  )»»M)»FACT(N)/(  (2.»''N)»FACT(N-H)»FACTIM1  1 

RETURN 
7  P=0.0 

RETURN 

END 

BOP   FACT  14 

FUNCTION  FACr(N) 
C      THIS  SUBPROGRAM  COMPUTES  H    FACTORIAL. 
FACr=l. 
IF(N) 1,2,4 
1  FACT=.999q99999999g99999999999999999999999999999999Et99 
RETURN 


^"1^ 


lOS 


<.    DO    5     1  =  1, N 

FACT=FACT«FLOAT( I ) 
5    CONTINUE 
2    RETUKN 

END 

BOP   BaLTZMANN2  1 't 

C      THIS  PROGRAM  IS  THE  CONTROL  PROGRAM  FOR  PHASE  2. 

DIMENSION  P0LY(6),EIGENV( 0 ) , H (  1 1 , 5) , GC ( 1 1 , 5 )  ,C ( 2 ) ,SI GMA ( 2 ) , E IGENC ( 
15),BC( Il.b) 

DIMENSION  B(  10), A( 10,10), IPC  10) 

COMMON  POLY,N,EIGENV,G,NREG,NBC,C,SlGMA,REFLTH,tlGENC,r,C,ACCUR,BC, 
1R,B,A,  IP 

REWIND  't 

REWIND  8 

00  35  L=l,6 

IP(L)=0 

B(L)=0.0 

DO  J5  K=l,6 
35  A(L.K)=0.0 

READ  TAPE  ^.NCASES 

DO  100  ICASE=1,NCASES 

WRITE  OUTPUT  TAPE  6,7,ICASE 

TYPE  7, ICASE 

7  FORMATI  16Xl'iHPR0BLEM  NUMBER,  15) 

READ  TAPE  'i ,  I  PROB  ,  POLY  ,  N  ,  E  IGFNV  ,  G  ,NREG  ,NBC  ,C,  S  IGMA  .REFLTH,  E  I  GENC.G 
IC,ACCUR,BC,R 

GO  TO  t8,22),IPRQB 
C      THIS  SECTION  SUPERVISES  SET  LP  OF  THE  PROPER  BOUNDARY  CONDITION 
C        MATRIX  IN  BC- 

8  GO  TO  (20,21,2't,  22)  ,NBC 

20  CALL  MARSHK 
GO  TO  22 

21  CALL  MARKBC 
GO  TU  22 

2A  IFIN-'t)  18,  18,  19 

ly  WRITE  OUTPUT  TAPE  6,17,N 

17  FORMAT ( IXbHTHE  P , I  3, 1X80HVAR I  AT  I ONAL  BOUNDARY  CONDITIONS  CANNOT  BE 
1  USED  AS  THEY  HAVE  NOT  BEEN  WORKED  OUT.) 

IPRUB=2 
GO  TO  22 

18  CALL  VARIHC(C(NREG) ) 

22  CALL  CHECKI 

100  WRITE  TAPE  8 , I PROB , POL Y , N ,E I GENV , G, NREG , NBC , C , S I GMA, REFLTH ,E IGENC, 
1GC,ACCUR,HC,R 

REWIND  'i 
REWIND  0 

WRITE  TAPE  A,NCASES 
DO  101  ICAS£=1,NCASES 

READ       TAPE    8, I PROB , POL Y , N ,E IGENV , G.NREG , NBC , C , S I GMA, RFFLTH , E IGENC , 
ICCACCUR.BC.R 

101  WRITE    TAPE    A , I PROB, POL Y , N , E I GENV , G, NREG , NBC , C , S I GMA, REFLTH , E I GENC , 


IGC,ACCUK,HC,K.8, A, IP 
REWINO    4 
CALL    EXIT 

srriH 

END 

BDP   MARSHK  I 'i 

SUBROUTINE  MARSHK 
C      THIS  SUBPROGRAM  SETS  UP  THE  COEFFICIENT  MATRIX  BC  WHICH  REPRESENTS 
C        MARSHAKS  OOUNOARY  CONDITIONS  AT  A  VACUUM  INTERFACE. 

DIMENSION  P0LYI6),EIGENVl5),n( 1 1 , 5) , GC ( 1 1 , 5)  ,C I  2 ) , SI GMA ( 2 ) .EIGENCI 
15),8C(  11,5) 
COMMON  PQLY,N,EIGENV,G,NREG,NBC,C,SIGMA,kEFLTH,EIGENC,r,C,ACCUR,BC 

N0V2=(Ntl)/2 
NP1=N+1 
DO  5  J=1,N0V2 
DO  b  LPl=l,NPl 
5  BCILPl, J)  =  BCMARSILP1-1,  J) 
RETURN 
END 

HOP   BCMARS  I't 

FUNCTION  BCMARSIN.J) 
C      THIS  SUBPROGRAM  COMPUTES  THE  VALUE  OF  THE  ANALYTICALLY  INTEGRATED 
C         INTEGRAL  FROM  -1  TO  0  OF  P(N,X)  TIMES  X  TU  THE  2J-1  POWER  FOR 
C        THE  MARSHAK  BOUNDARY  CONDITIONS  AT  A  VACUUM  INTERFACE. 

IFINJI, 1,2 

1  BCMARS=-1./FL0AT(2«J) 
RETURN 

2  BCMARS^O.O 
MP1  =  N/2H 

DO  b  KPl=^l,MPl 
K=KP1-1 
5  BCMARS=BCMARSt((-l.)»»K»FACTI2«N-2«K) I / ( 2. ••N»FAC T ( K ) »FACT ( N-K ) "F A 
ICTIN-2«K1«FLQAT(N+2»J-2»K)) 
BCMARS=BCMARS«(-l.)»»INtl) 
RETURN 
END 

BOP   MARKBC  14 

SUBROUTINE  MARKBC 
C      THIS  SUBPROGRAM  SETS  UP  THE  COEFFICIENT  MATRIX  BC  WHICH  REPRESENTS 
C         MARKS  BOUNDARY  CONDITIONS  "T  A  VACUUM  INTERFACE. 

DIMENSION  P0LY(6),EIGENV(  !j),G(11,5),GCI11,5),C(2)  ,S1GMA12)  ,EIGENC( 
151,HC( ll.b) 

COMMON  P0LY,N,EIGENV,G,NREG,N8C,C,SIGMA,REFLTH,E1GENC,GC,ACCUR,BC 

N0V2=lNtl)/2 

NP1=N« I 

K=l 

X=-1.0 

ANS1=PNP1(X) 

DO  !)  1  =  1,39 
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X=Xt0.025 
ANS2=flNSl 
ANSl  =  PiMPl(X) 

IFISIGNK  1., ANSI) -SIGN  I l.,ANS7) ) If  5,1 
F      PNPl 

1  VALUE=ROOT(PNP1,X-0.02  5,X,0.,ACCUR) 
DC  2  J=1,UP1 

2  BC( J,K)=P(N-NPltJ,VALUE) 
K  =  K  +  1 
IF(K-N0V2I5,5,3 

'j  CCNriNUE 

3  KETURiN 
END 

BOP   PNPl  I'. 

FUNCTION   PNPl(X) 

DIMENSION  P0LY(6) 

COMMON  POLY.N 
C      THIS  SUBPROGRAM  IS  USED  BY  THE  MARKBC  SUBPROGRAM  TO  FIND  THE  ROOTS 
C        OF  THE  EQUATION  P(N+l,X)=0. 

PNPl=P(Ntl,X) 

RETURN 

END 

nop   VARIBC  I'l 

SUBROUTINE  VARIBCICP) 
C      THIS  SUBPROGRAM  SUPERVISES  THE  SET  UP  OF  THE  BOUNDARY  CONDITION 
C        MATRIX  8C  WHICH  REPRESENTS  THE  VARIATIONAL  BOUNDARY  CONDITIONS. 
DIMENSION  P0LY(6),EIGENV(5),n( 1 1 , i I , GC t 1 1 , 5 )  ,C ( 2 ) ,SIGMA(2) ,EIGENC( 
15),BC(  n.-j) 
DIMENSION  COC,'.) 

COMMON  P0LY,N,EIGENV,G,NREG,NBC,C,S1GMA,REFLTH,EIGENC,GC,ACCUR,BC, 
IR.CCCPRIME 
IF(N-2)  1,2,'. 
I  BC(  l,l)=S(jRT(2.)/3. 
GO  TO  3 


UP  THE  NECFSSARY  COEFFICIENTS  TO  FIND  THE 


2 

RC(  1,1)  =  SURT(2.)/ 

i 

BC(2,1)=-1. 0/3.0 

RETURN 

C 

THIS  SECTION  SETS 

c 

P-3  CONSTANTS. 

'i 

BCl^, 11=1.0/7.0 
BC(  3,21  =  1.0/5.0 
IFlN-3)5,5,6 

5 

CU(  l,l)  =  Ui. 
CQI  1,2)=-A0. 
CC(  1,31  =  100. 
C0(  1,'>1=-117. 
CQ(2, 11=20. 
C0(2,21=-70. 
CC(2,3)=',9. 
CC(2,^1=-16. 
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CU( 3,11=1. 

C0( 3,2)=-^. 

C0( 3,3)  =  ). 

CC(  i,',)'-.U 

CG(^,1 )=-72. 

CC('.,?)  =  S6. 

CG(',,3)  =  100. 

CCl'.,'i)=-l'.7. 

CPfUMC  =  0.0 

GC  TU  7 
C      THIS  SCCriON  SETS  UP  THE  NECESSARY  COEFFICIENTS  TO  FIND  THE 
C        P-4  CONSTANTS. 

6  CQ(l,l)  =  16.«(9.  +  (  l.-CPl»(2.+<5.«(  l.-CP)  )  ) 
GUI  l,2)=-2.«(  l'.0.-4S0.»(  l.-CP)  ) 

CO  (  1,3  1  =  2'.  50. 

CC( 1,^)=-1A7.»9. 

CQ(2,l)  =  161.«(l't0.-'t90.*(l.-CPl)  +  32.*2450.»(l.-CP) 

CGI 2,2)=-2.«9.#1715. 

CG(2,3)=27.»147.«9.+32.*2  52.»( l.-CP) 

C0(2,',)=-161.»l'.4.-301.»32.»(  l.-CP)-469.»16.»(  l.-CP)  ••? 

CGI  3, l)  =  27. 

CGI 3,2)=-2.«27. 

CG( 3,3)=161. 

C0(  3,',)=-161. 

CPR1M£=32.«( l.-CP) 

C0(4,  l)=-6'.8. 

C0(4,2)=7.»9.«8. 

C0(  4,3I=2'.50. 

CG(4,',)  =  -IA7.»9. 
C      THIS  SECTION  SUPERVISES  THE  FINDING  OF  THE  CONSTANTS  A, B.C.  AND  D 
C        FOR  THE  P-3  AND  P-'t  APPROXIMATIONS. 

7  E  =  .5 
DC=.l 

10  ANSl=FQURfca(E) 

8  ANS2=ANS1 
IFIBC(2,  1)  113,14,1', 

13  ITEST=1 
GO  TO  15 

14  ITEST=0 
E=E+DC 

ANSl=FOURCQ(e) 

IF(SIGN(  l.,ANSl)-SIGNI l.,ANS2) ) 9,  11,9 
IF(E-2. 0)8, 0,12 

12  E=.5 

DC=DC/10. 
GG  TO  10 
FGUREG 

9  IF( ITEST)H,8, 16 
16  F=R0OT(FQUREQ,E-DC,E,O.0,ACCUR) 

CALL  NQNLINIF) 
IF(liC12,  1)  )  17,  10,  18 
18  CALL  NQNLINI E) 


15 


1  I 
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GO  TO  8 
17  BCI2,ll  =  BC(2,n/3.0 
BC(2,2)=RC(2,2)/3.0 

RETURN 
END 

BOP   FOUREQ  I'l 

FUNCTION  FOUREQ(E) 
C      THIS  SUBPROGRAM  IS  USED  AS  AN  AUXILLIARY  PROGRAM  BY  VARIBC  IN  THE 
C        FINDING  IF  THE  CONSTANTS  A.B.C,  AND  D  FOR  THE  P-3  ANO  P-4 
C        APPROXIMATIONS 

CALL  NONLINIE) 

FOUREQ=E 

RETURN 

END 

BOP   NONE  IN  I'l 

SUBROUTINE  NONLIN(E) 
C      THIS  SUBPROGRAM  SOLVES  THREE  OF  THE  FOUR  NONLINEAR  EQUATIONS 
C        AND  USES  THE  FOURTH  AS  A  CPNSISTENCY  CHECK. 

C      IN  THE  USUAL  NOTATION  A  =  BC ( 1 , 1  I , B=BC ( 2 , 1 1 ,C  =  BC ( 1 , 2) , AND  D=BC(2,2). 
DIMENSION  P0LY(6),EIGENV(5),'-(  II  ,  5  )  ,  GC  I  1 1 ,  5  )  ,C  (  2  )  ,  SIGMA!  2  )  ,EIGENC( 
15),BC(U,'J) 
DIMENSION  Ca(^,4) 
COMMON  P0LY,N,EIGENV,G,NRCG,N8C,C,SIGMA,REFLTH,EIGENC,GC,ACCURfBC, 

IR.CO.CPRIME 
BC( 1.2)=E 
BC(  1,1)  =  S0RT(  IC0(1.1)  +  BC1  l,2)»{COll,2)+BC(l,2)»CO(l,3)))/(-C0(l,'.) 

n ) 

RCI2,2)=-IBC(l,n»(C0l2,3)+Br( 1,2)«C0(2,2)) )/(CQ(2,l)»nc(l,2)+C0(2 

l.A)  ) 
BC(2,l)='-(C0(3,n+BCl  l.2)»C0(3,2)*CU(3,'.)«BCll,n»BC(?.2))/(BC(l.2 

l)«C0(3,3)tCPRIM£) 
E  =  C0('.,l)  +  BC(2,l  l»(C01'i,2)+BC(2,l)»C01't,'t)  )+BC(2,2)»0CI2,2)»C0l4,3 

I) 

RETURN 
END 

OOP   P  1^ 

FUNCTION  P(N,X) 
C      LEGENDRE  POLYNOMIAL  FUNCTION  CALCULATOR   , 
IF(N)1,2,J 
100  FORMATI 1HB,22HNEGATIVE  ORDER  N  FOR  P) 

1  WRITE  OUTPUT  TAPE6,100 

2  P=l. 
RETURN 

3  M=  N/2 
IFIXj^.-itb 

4  IF(N-2»M)7,8,7 

a  P=(  (-1. )««M) 'FACTIN)/! (2.»«N)»FACT(N-M)«FACT(M)  ) 

RETURN 
7  P=0.0 
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RETURN 

5  IF(N-1 )2, II,  16 

16  TERM=(  (-1.  1»»M)»FACT(2»N-2*M)»(X'>«(N-2»M)  )/(  (  2  .  ""N) 'F  ACT  (  M)  »FACT  (  N 
1-M)  ) 

SUM=TERM 

i;)06K=l,M 

TERM=-TERM«FL0AT(M-K+l)»FLOAT(2»(N-M+K)-U»FLUAT(2»(N-M+K) )«(X»»2) 
l/(FL0Ar(N-K+K)«FL0AT(N-2»IM-K)-l ) "FLOAT ( N-2» I M-K ) ) ) 

6  SUM^SUM+TERM 
P=SUM 
RETURN 

11  P  =  X 

RETURN 
END 

HOP   BQLT2MANN3  lA 

C      THIS  PROGRAM  IS  THE  CONTROL  "ROGRAM  FOR  PHASE  A 

DIMENSION  PULY(6),EIGENV(  'i  )  ,  C  (  1 1  ,  i  )  ,  GC  (  1 1 ,  5)  ,C  I  2  )  .SIGMA  (2)  .EIGENCI 
15),BC(ll,b),A(10,101,IP(10),V(lO).R(10) 

COMMON  P0LY,N,E1GENV,G,NREG,NBC,C,SIGMA,REFLTH,CIGENC,GC,ACCUR,BC, 
IR,A,IP,V,B,NSUB 
REWIND  't 
REWIND  8 

READ  TAPE  A,  NCASES 
ITIWE=1 

DO  21  KCASE=1, NCASES 
WRITE  OUTPUT  TAPE  6,3,KCASE 
TYPE  3,KCASE 
3  FCRMAT( 16X14HPRQBLEM  NUMBER, 15) 
KrAPE=4 

1  READ  TAPE  KT APE,  I  PROB , POL Y,N . E IGENV.G ,NREG , NBC ,C , SIGMA .REFLTH, E  I  GE 
INCGC,  ACCUR,BC,R,B,A,  IP 

GC  rO  (C,22) , ITIME 
8  GO  TO  (55,20) , IPROB 

2  WRITETAPE  KTAPE, I PRUB, POL Y,N , E IGENV i G , NREG.NBC ,C , SIGMA .REFLTH, E I GE 
1NC.GC,ACCUR,BC.R.B,A, IP 

GC  TO  121,23) .ITIME 
C      THIS  SECTION  SUPERVISES  THE  SETTING  OF  THE  CRITICALITY 
C        DETERMINANT  TO  ZERO  THUS  FINDING  A  CRITICAL  RADIUS. 
55  0R=.2«R 

ANS1=CRITEQ(R) 

15  ANS2=ANS1 
R=R*DR 

ANS 1=CRITEQ(R) 
IFI ANS1«ANS2) 16, 16, 15 
F      CRITEO 

16  R=ROOT(CRITEO,R-DR,K,0. ,ACCUR) 
CALL  SETUPA(R) 

20  KTAPC  =a 
GO  TO  2 

21  CONTINUE 
REWIND  ', 
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REWINO    8 

WRITE    TAPE    A.NCASES 

ITIME=2 

DC  23  KCASE=l,NCASeS 

KTAPE=8 

GO  \0    1 

22  KTAPE='. 
GO  TO  2 

23  CONTINUE 
REWIND  A 
CALL  EXIT 
STOP 

END 

HOP   CRITEQ  I'l 

FUNCTION  CRITEQ( R) 
C      THIS  SUBPROGRAM  IS  USED  BY  B0LTZMANN3  AS  AN  AUXILLIARY  PROGRAM 
C        TO  HELP  IN  FINDING  THE  CRITICAL  RADIUS, 

DIMENSION  P0LY(6),EI6ENVC>),n(  U  ,  5  )  ,  GC  (  1 1 ,  5  )  ,C  (  2  )  ,  SI  GMA  (  2  1  ,  E  IGENC  ( 
15),DC{11,^),A(10,10),IP(10),V(10),B(10) 

COMMON  POLY,N,EIGENV,G,NREG,NBC,C,SIGMA,REFLTH,EIGENC,GC,ACCUR,BCi 
1Z,A, IP.V.B.NSUB 
CALL  SETUPAIR) 
N0V2=(N+l)/2 
IFINBC-'i)  16,  15,16 

15  NSUB=2«N0V2 
GO  TO  17 

16  NSUB=(U2»(NREG-1)  l«N0V2  , 

17  CALL  CRAM(NSUB,1) 
CRITEO=DET(NSUa) 
RETURN 

END 

BOP   SETUPA  I'l 

SUBROUTINE  SETUPA(R) 
C      THIS  SUBPROGRAM  SETS  UP  THE  CRITICALITY  DETERMINANT  IN  A. 

DIMENSION  P0LY(6),EIGENV(5),n( ll,5),GClll,5)  , D ( 2 ) , SI GMA ( 2  I ,E IGENC ( 
I5I,BC( 11,5),A( 10, 10) 

COMMON  POLY,N,£IGENV,G,NREG,NBC,D,SIGMA,REFLTH,EIGENC,GC,ACCUR,BC, 
IRA, A 
DC  24  1  =  1,  10 
00  24  K==l,  10 

24  A( I  ,K)=0.0 
NP1=N+1 
N0V2=NPl/2 
N0V=2»N0V2 

DO  10  K=l,N0V2 
ARG=SIGMA( l)»R/EIGeNC(K) 
IF(K-l)  I,  1.2 

1  1  =  1 

GO  ro  3 

2  1  =  0 
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J  TEMPI=C(0,ARG, I) •GC( 1,K) 
DO  10  LP1=1,NP1 

TEMP2  =  C(LPl-l,AKG,I)»GC(LPl,K)»FLaAT(2»LPl-n 
GG  TO  i4,6),NREG 
C      THIS  SECTION  SETS  UP  THE  MATRIX  FOR  A  BARE  SPHERICAL  REACTOR, 
'i  UO  b  NSUB=1,N0V2 

b  A{NSUB,K)=A(NSUB,K)+BCILPI,NSUB1»TEMP2 
CC  TO  10 
C      THIS  SECTION  SETS  UP  THE  PARTS  OF  THE  MATRIX  KCPRESCNTING  THE 
C        INTERFACIAL  BOUNDARY. 

6  IF(N-2»(N/2) )7,8,7 

7  A(LPl,K)=rEMP2 
GO  TO  10 

8  IF(LPl-l) 10, 10,9 

y  A(LP1-1,KI=TEMP2-P(LP1-1,0.0)»TEMPI»FLOAT(2»LP1-1) 

10  CONTINUE 

GO  TO  (19,11) ,NREG 

11  DO  16  1=1,2 

DO  16  K=1,N0V2 

ARG= (-!.)••( I-1)»SIGMA(2)»R/FIGENV(K) 
ir(igBC-4)22,20,22 
2U  IF( I-l) 16,16,21 

21  NSUB=K+NQV2 
GO  TO  23 

22  NSUH=KtN0V2»I 

21  TEMP1=QI0, ARG)»G( 1,K) 
DC  l'>  LP1=1,NP1 

TCMP2=-0(LPl-l,ARG)»FLOAT(2»LPl-l )«G ( LPl , K ) • (-1. ) ••( I LPl-l ) • ( I-l ) ) 
IF(N-2«(N/2) ) 12, 13,12 

12  AILP1,NSUB)=TEMP2 
GO  TO  15 

13  IF(LP1-1I 15, 15,14 

14  A{LPl-l,NSUB)=TEMP2+P(LPl-l,0.0)»TEMPl»FLaAT(2»LPl-l) 

15  CONTINUE 

16  CONTINUE 
1FINUC-4)17, 19,19 

C      FOR  A  REFLECTED  REACTOR  THIS  SECTION  SETS  UP  THE  VACUUM 
C        INTERFACE  BOUNDARY  CONDITIONS. 

17  DO  18  1=1,2 

DO  IB  K=1,N0V2 

ARG=(-1. )••(  I-1)«SIGMA(2)»(R+REFLTH)/EIGENV(K) 

DO  18  LP1=1,NP1 

TEMP2=Q(LPl-l,ARG)»FLOAT(2»LPl-l)*G(LPl,K)»l-l.)»»(ILPl-ll»(I-n) 

NSUB  =  N0V2»  It-K 

DO  IB  L=1,N0V2 

NSULU  =  L  +  NUV 

18  A(NSUtU,NSUB)  =  A(NSUBl,NSUB)tRC(LPl,Ll»TCMP2 

19  RETURN 
END 

BOP   CRAM  14 

SUBRUUTINE  CRAM( N, I ) 


')\r 
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C      CROUT  REDUCTION  OF  AUGMENTED  MATRICES 

C      THIS  PROGKAM  PERFORMS  A  CROUT  REDUCTION  ON  A  MATRIX  A. 

C      WITH  1=1,  THE  CROUT  REDUCTION  IS  PERFORMED  WITH  ROW  INTERCHANGES. 

C      WITH  1=2,  THE  CROUT  REDUCTION  IS  PERFORMED  WITHOUT  ROW  CHANGES. 

DIMENSION  P0LY(6),EIGENVI5),G( I 1 , 5) , GC ( U . 5 ) ,C ( 2 1 , SIGMAI 2 ) .EIGENCl 
15),BC(11,5),A(10,10),IP(101,VI10),BI10) 

COMMON  POLY,M,EIGENV,G,NREG,NBC,C,SIGMA,REFLTH,EIGENC,GCtACCUR,BC, 
1R,A, IP,V,8,NSUB 
GO  TO  (2200,2201  )  , I 

2200  IDMV=1 
G0TQ22O2 

2201  IDMV=2 

C      REDUCTION  OF  MATRIX 

2202  IF(N-l)2223,2223,222't 
2223  IP( l)=0 

RETURN 
222'>    DO  2204  IOK=l,N 

V( IDK)=ABS(A( IDK, 1) ) 

D0220'.  101=2, N 

IFI VI I0K)-ABS(A( IDK, IDI ) ) ) 2203, 2  204, 2204 

2203  V(IOK)=ABS(A(IDK,IDI)) 

2204  CONTINUE 

DO  2222  IUK=1,N 

DETR=-1. 

I0K1=IDK-1 

DU2214  IDI=IDK,N 

DCTPK=0.0 

IFI IDK- I) 2208, 2208, 2206 

2206  DU2207  IDJ  =  1, I  OK  1 

2207  DETPR=DETPRtA(IDI,IDJ)«A(IDJ,IOK) 

2208  A{ IDI, IDK)=A( IDI , I0K)-OETPR 
GQ  Ta(22l2,2225) , IDMV 

2212  DETS=ABS(A(IDI,IOK) )/V( IDI) 
IFI DETS-DETR)2214, 2214,  2213 

2213  DETR=DETS 

IP(  IDK)  =  IIJK-IOI 
GO  TO  2214 
222i  IP( IOK)=0 

2214  CONTINUE 
IDK2=1DK-1P( IDK) 
DETR=A( IDK2, IDK)/V( IDK2) 

IFI ADS (DETR)-l.E-08 12230,2230,2232 

2240  FORMAT!  lHB,t)HPlVQT,  13,  ISHIS  LESS  THAN  l.E-08) 

2241  FOKMATI 1HB,^HP1V0T, I3,7HIS  ZERO) 
2238  WRITE  OUTPUT  TAPE  6,2240, IDK 

IFI A(  IDK2,  IDK) )  2232,2231,2232 

2231  WRITE  OUTPUT  TAPE  6,2241, IDK 
IF  I N-NSUB 12501, 2500, 2500 

2501  CALL  EXIT 

2232  V( I0K2)=VI IDK) 
VIIOK)=DETR 
DQ2222  I0J=1,N 
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DETR=A( IDK, lOJ) 

IF  (  IDJ-I  010  2215,  2215,2216 

2215  A( lOK, IDJ)=A( IDK2, lOJ) 
GO    ru    2220 

2216  DETPR=0.0 

IF(  IDK-1)2219,22I'3.2217 

2217  DQ2218  IDI  =  l,  I  L'K  1 

2218  DETPR=OET|JRtA(IDK,IDI)«AIIOI,IDJ) 

2210  AdDK,  IDJ)  =  (A(  IDK2,  I  DJ  )-OETPR  I /A  (  lOK.IOKI 

2220  IF(  IP(  lOK)  )2221, 2222,2222 

2221  A( I0K2, IDJ )=OETK 

2222  CCNTlNUe 
2500  RETURN 

END 

BOP   DET  14 

FUNCTION  UET(N) 
C      AFTER  CALLING  THE  CRAM  SUBROl'TINC  THIS  FUNCTION  WILL  COMPUTE 
C        THE  DETERMINANT  OF  THE  MATRIX  A. 

D I  MENS  I  UN  PULY(6),ElGENV(5),C(ll,5),GC(ll,5» ,C ( 2 » ,SIGMA(2 » ,E IGENCI 
L5),BC(  U,5)  ,A(  10,  10),  1P(  10)  ,V(  10)  ,H(10) 

COMMON  POLY,M,EIGENV,G,NREG,NBC,C,SIGMA,ftEFLTH,ElGENCtGCtACCURiBC, 
1R,A, IP,V,B 
0ET=1.0 

DO  2229  IDK=1,N 
DET=DET»A( IDK, lOK) 
IF(  IP(  IDK)  )2223, 2229,2229 

2223  DET=-DET 
2229  CONTINUE 

RETURN 
END 

BOP   C  lA 

FUNCTION  C(L,X, I  ) 
C      THIS  SUBPROGRAM  COMPUTES  THE  VALUt  OF  THE  FUNCTION  CIL.X)  FOR  BOTH 
C        REAL  AND  NEGATIVE  IMAGINARY  ARGUMENTS.  L  I S  THE  ORDER  OF  THE 
C        FUNCTION  AND  X  IS  THE  ARGUMENT  OF  THE  FUNCTION.  FOR  REAL 
C        ARGUMENfS  1=0,  FOR  NEGATIVE  IMAGINARY  ARGUMENTS  1=1. 
IFIX) 10,9, 10 
9  IFIL) 12,  12, II 

11  C=0.0 
RETURN 

12  C=l. 
RETURN 

10  IFI  D't.'t,  1 
4  C=.5»lQ(L,X)+(-l. )»*L»Q(L,-X)) 
RETURN 

1  IF(L-2»(L/2) )2,3,2 

2  COFUNC=-CUS(X) 
FUNC=-SIN(X) 
SIGNB^l. 

GO  10  7 
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COFUNC=SINtX) 

FUNC=CQSIX) 

SIGNn  =  -l. 

C=(CUrUNC.CSER(L,X,0)-SlGNR»FUNC»CSER(L,X,l) )/X 

RETURN 

END 


BOP   CSER  M 

FUNCTION  CSER(L,X,  I  ) 

THIS  SUBPROGRAM  COMPUTES  THE  FUNCTIONAL  VALUE  OF  THE  SERIES 

NCCESSAKY  FOR  THE  FUNCTION  C(L,X).  THE  ARGUMENTS  OF  L  AND  X  ARE 
THE  SAME  AS  IN  THE  FUNCTION  CIL.X.l)  SUBPROGRAM.  1=0  FOR  A 
STARTING  INTEGER  OF  0,  AND  1=1  FOR  A  STARTING  INTEGER  OF  I. 

CSER=0.0 

IF(L-I )3, 10, 10 
10  TERM=FACT(L+I )/( FACTl Il*FACr(L-n»(2.«X)»»I) 

CSCR=TERM 

IP2=1+Z 

IF(L-IP2) 3,1,1 

1  rMULT  =  -.2'j/(X»Xl 

DO  2  J=IP2,L,2 

TERM  =  TERM«FL0AT(L-J<-2)«FL0AT(L-JtU»FL0AT(L  +  J-Ll*FL0AT(LtJl»FMULT/ 

I (FLOAT( J-1 ).FLOAT( Jl ) 

2  CSER=CSER+TERM 

3  RETURN 
END 


e.OP   Q  I'l 

FUNCTION  fJlL.X) 

THIS  SUBPROGRAM  COMPUTES  THE  VALUE  OF  THE  FUNCTION  Q1L,X) 
Q=(CXP(X1 )/X 
I  F  (  X  )  7  ,  7  ,  R 

7  SICNA=l. 
X--X 

GO  TO  9 

8  SIGNA=-1. 

9  P0LY=1. 
FMULT=.5»SIGNA/X 
AN=L+1 
ANPLUS^L 
TERM^l. 

N=l 
k    TERM=TCRM»FMULT«ANPLUS*AN/FLaAT(N) 

POLY=POLYtTtRM 

IF(  L-NI6,6,'i 
ii  ANPLUS  =  ANPLUS-l. 

N  =  Ntl 

AK=AN<-  1. 

CG  TO  -i  • 

6  a=0»POLY 
2  RETURN 

END 
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POP   8QLTZMANN'.         14 
C      THIS  PRUGRflM  IS  THE  CONTROL  PROGRAM  FOR  PHASE  A 

DIMCNSION    PULYIfij.CIGENVCjl.GI  U,5),GCIU,5)  ,C(2),SIGMA(2)  ,EIGENC( 
L5),HC(11,5),A(10,10),IP(10),V(10),B(10),S(10,10),T(10) 

CCMMUN  PULY,N,EIGCNV,G,NREG,NBC,C,SIGMA,REFLTH,EIGENCtGC,ACCUR,BC, 
IRiA,  IP,V,B,l\iSUB,  S,T 
PI=3. 14159  265  3  58979324 
REWIND  4 
REWIND  8 

READ  TAPE  4,  NCASES 
DO  1000  ICASC=1, NCASES 
WRITE  OUTPUT  TAPE  6,7,  ICASE 
TYPE  7,  ICASE 
7  FORMAT! 16X14HPRUBLEM  NUMBER, 15) 

READ  TAPE  4, I PROR , POLY, N, E IGFNV , G ,NREG ,NBC ,C , S IGMA .REFLTH, E I GENC ,G 
IC,ACCUR,BC,R,B,A,  IP 
GU  TO  (8, 35) , IPRQB 
C      THIS  SECTION  ASSUMES  Al  IS  ONE  AND  COMPUTES  THE  VALUES  OF  ALL 
C        THE  OTHER  CONS TANTS , A-SUB  K. 
B  N0V2=(N*l)/2 

IF(NBC-4) 10,9, 10 
9  NSUB=2»N0VZ 
GO  TO  11 

10  NSUB=( 1+2«(NREG-1 ) )»N0V2 

11  iriNSUB-l)  12,32,  12 

12  NSUH1=NSUB-I 
NSUB2=NOV2»NREG 
DO  IB  L=1,NSUB1 
1F(L-NSUB2) 17, 15, 15 

15  00  16  K=1,NSU8 

16  A(L,K)=A(L+1,K) 

17  B(L)=^-A1L,  1) 

DO  18  K=l,NSUBl 

18  A(L,K)=A(L,K+1) 
DO  19  K=l,NSUBl 
T(K)=B(K) 

DO  19  L=1,NSUB1 

19  SIL,K)=AiL,K) 
CALL  CRAMINSUBI, 1) 
CALL  SOLVE (NSUBl ) 
CALL  RESIDUINSUBl  ) 

c     THIS  secriCN  normalizes  the  total  neutron  flux  so  that  AT  R  =  0 
c      IT  is  one  neutron  per  souare  centimeter  per  second. 

31  DO  30  K=l, NSUBl 
L=NSUBlt2-K 

30  B(L)=B(L-1) 

32  B(l )=1.0 
SUM=0.0 

DD  33  K=1,NUV2 

33  SUM=^SUM  +  B(K) 
SUM=SUM»4.«PI 


lis 

"^-    .... 


SUM=SUM-5.•ACCUK•SUM 
DQ    i'f    K=1,NSUB 
i'l    B1K)  =  B(K)/SUM 


35    CALL    CHECKI  cocr-    r 

1000    WRITETAPE    8, I PKOB , POLY , N, E IGENV, G ,NKEG ,NBC ,C , S JGMA .REFLTH, EI GENC , b 

ICtACCUK.BCR.B 


ICtACCUR.BCR.B 
REWIND    A 
REWIND    8 

WRITE  TAPE  4,  NCASES 

DO  36  ICASE  =  1, NCASES  ^,,  ^,^c.,r  r 

READ  TAPE  8, I PROB , POLY , N, EIGENV , G ,NREG ,NBC .C , S IGMA .REFLTH, E IGENC ,G 

1C,ACCUR,BC,R,B  ^,  ^  ,  ^.^^„^  ^ 

36  WRITETAPE  A, I PROB , POLY, N, E IGFNV, G ,NREG , NBC ,C , S IGMA .KEFLTH.E IGENC ,G 

lC,ACCURtBC,R,B 
REWIND  't 
CALL  EXIT 
STOP 
END 

BOP   SOLVE  lA 

SUBROUTINE  SOLVE  (N)  ^ , ^^ 

C      AFTER  CALLING  CRAM  THIS  SUBROUTINE  WILL  COMPUTE  THE  SOLUTION 
C        VECTOR  UF  THE  MATRIX  EQUATION  AX=B.  BEFORE  RETURNING  TO  THE 
C        MAIN  PROGRAM  THE  SOLUTION  VECTOR  IS  STORED  IN  B. 

DIMENSION  P0LY(6),ElGENV(5),n(  ll.5),GCm,5)  ,C  1  2  )  ,  SI  GMA(2  I  .EIGENC  ( 
15),BC(ll,bl,A(10,10),IP(lO).V(lO),B(lOI 
COMMON  POLY,M,ElGENV,G,NREG,NBC,CfSIGMA,REFLTH,CIGENCtGCiACCURiBCi 

IR.A, IP.V.B 
DQ  2256  IDK  =  1,N 
IDKl  =  IDK  -  1 
IDK2  =  IDK-IP( IDK) 
DETR  =  Bl IDK) 
DETPR  =  0.0 

IF( IDK-l)  2253,  2253,  2257 
2257  DQ  2252  IDI  =  I,  I DK I 

2252  DETPR  =  DtTPR+A( I DK, I D I ) •  B(IOI) 

2253  B(IDK)  =  (B(IDK2)  -  DETPR)  /  AIIOK.IDK) 
IF  IIP(IDK))  2254,  2256,  2256 

225'>  B(IDK2)  =  DETR 
2256  CONTINUE 

DO  2263  IDI2  =  l,N 

lOI  =  N  +  1  -  I0I2 

DETPR  =    0.0 

IDI  1=101  +  1 

IF  (N  -  IDI)  2263,  2263,  2261 

2261  DO  2262  ICJ  =  IDIl,N 

2262  DETPR  =  DETPR  +  A(IDI,IDJ)»  BIIOJ) 

2263  B(  IDI  )  =  BdDI  )  -  DETPR 
RETURN 

END 

BOP   RESIDU  14 


SUBKOUTING  RESIDU(N) 
C      AFTER  THE  CRAM  AND  SOLVE  SUBROUTINES  HAVE  BEEN  CALLED  THIS 
C        SUBROUTINE  WILL  COMPUTE  TH»=  RESIDUALS  IN  THE  COEFFICIENT 
C        VECTOR  T  AND  ITERATE  ON  THE  ANSWER  VECTOR  D  UNTIL  THFRE  IS 
C        NO  CHANGE  IN  8  FROM  ONE  ITFRATION  TO  THE  NEXT.  THE  SUBROUTINE 
C        ASSUMES  THAT  THE  ORIGINAL  t'ATRIX  IS  IN  S  ANU  THAT  THF  ORIGINAL 
C        COEFFICIENT  VECTOR  IS  IN  T. 

DIMENSION  P0LY(6),EIGENV(5),C(  ll,t>),GC(ll,5),C(2)  ,SIGKA(2)  tEIGENCI 
15).BC(U,5),A(10,lO)fIP(10),V(lO).BI10),S(10,ia),T(10) 

COMMON  PQLY,M,EIGENV,G,NREG,NBCiC,SIGMA,REFLTH,CIGENCtGC,ACCUR,BC, 
1R,A,1P,V,B,NSUB,S,T 

DO  1  1=1, N 

VII  )«B(  I) 

DO  1  J=1,N 

1  T(I  )  =  T(  I)-S(I,J)<B(J) 

2  DO  3  1=1,N 

3  B( I)=T( I ) 
CALL  SOLVE  (N) 
DO  10  1=1, N 

DO  10  J=1,N 
10  T(I )=TI I)-S( I,J)»B( J) 
J  =  0 

DO  5  1=1, N 
B(I )=B( I )tVI I) 
IFIBd  )-V(  I  )  )6,5,6 

6  J  =  l 

5  V(I )=B(  I) 
IF(J)  7,7,2 

7  RETURN 
END 

OOP   B0LTZMANN5  l-i 

C      THIS  PROGRAM  ISTHE  CONTROL  PROGRAM  FOR  PHASE  5 

DIMENSION  POLY (6) , E I GENV ( 5 ) , n ( 1 1 , 5 ) , GC ( 1 1 , 5)  ,C ( 2 1 , SI GMA (2 ) ,E IGENC ( 
15),BC(  ll,ti),B(  10) 

DIMENSION  X( 100) ,Yl 100) 

COMMON  POLY,N,EIGENV,G,NREG,NBC,C,SIGMA,REFLTH,EIGENC,r,CfACCUR,BC, 
1R,B, IPROB 

COMMON  X,Y 

REWIND  4 

REWIND  8 

READ  TAPE  'i.NCASES 

DO  16  ICASE=1,NCASES 

WRITE  OUTPUT  TAPE  6,7,ICASE 

TYPE  7,  ICASE 
7  FORMAT!  16X1'.HPR08LEM  NUMBER,  T5) 

READ   TAPE  4, I PROB, POLY ,N , E I HENV , G,NREG , NBC , C , S IGMA , REFLTH , E IGENC , 
1GC,ACCUR,BC,R,B 

14  CALL  INTOUI 

15  CALL  CHECKI 

16  CONTINUE 
REWIND  4 
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REWIND    8 

WRITE    TAPE    A.NCASES 

DO    17    ICASE=l,NCASES 

READ   TAPE  8,  I  PROS, POLY, N ,E I CENV , G.NREG , NBC , C , S 1 GMA, RCFLTH , E IGENC. 

IGCACCUR.HC.R.B  ,-    ^^    ^ 

WRITE    TAPE    4, I PROB, POLY ,N ,E IRENV , G.NREG , NBC , C , SIGMA, REFLTH, EIGENC , 

1GC,ACCUR,BC,R,B 

GO  TO  (  18,17) , IPROB 
18  READ  TAPE  8,X,Y 

WRITE  TAPE  '(,X,Y 

READ  TAPE  8,iMSUB 

WRITE  TAPE  A,NSOB 

IF(NSUB)17,17,20 
20  DO  33  K=l,NSUB 

READ  TAPE  8, XM IN , XMAX, YMIN, Y^AX , I  SCALE , JSCALE , NPTS.NPLOTS , NCOPY, NC 

WRITE  TAPS',,  XMIM,XMAX,YMIN,Y^'AX,  I  SCALE,  JSCALE,NPrS,NPLOTS,NCOPY,NC 

lARDS 

READ  TAPE  8,X,Y 
33  WRITE  TAPE  'i,X,Y 
17  CONTINUE 

REMIND  4 

CALL  EXIT 

STOP 

END 

BOP   INTUUT  I'l 

SUBROUTINE  INTOUT 
C      THIS  SUBROUTINE  COMPUTES  ALL  OF  THE  NECESSARY  DATA  FOR  THE 
C        PRINTOUT  OF  THE  FINAL  RESULTS.  .,^^.,^, 

DIMENSION  P0LY(6),EIGENV('>),n(  1 1  ,  b  )  ,GC  (  1 1 ,  5  )  ,C  (  2  )  ,  SIGMA  (  2  )  ,EIGENC( 

15),BC(  ll,!il,fll  10) 

DIMENSION  X(100),Y1 100)  

COMMON  POLYiN,EIGENV,G,NREG,NBC,C,SIGMA,REFLTH,EIGENC,GC,ACCUR,BC, 

1R,B, IPROH 
COMMON  X,Y 
GO  TO  I  I't,  13)  ,  IPROB 
C      THIS  SECTION  COMPUTES  THE  DATA  FOR  A  NORMAL  PRINT-OUT 
I'l  NP1  =  N+1 
NMAX^lOO 
DO  101  K=l,NMAX 
101  Y(K)=0.0 

PI  =  3. 14159  265358  97932'! 

ISCALE^O 

JSCALE'O 

NC0PY=1 

NCARDS=3 

IFfNBC-'.ISS, 54,55 

54  NREG=l 

55  K0UTER=R*FL0AT(NREG-1)»REFLTH 

DX=-.l 
X(l)=1.0 
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19 

56 

13 

1 

15 


DC  102  K=2,21 

X(K)=X(K-1)+DX 

DO  19  L=1,NP1 

TEMP=PHIL(L-1,R0UTER,NREG) 

DO  19  K=l,21 

Y(K)=Y(K)+P(L-l,XtK))»TEMP 

PCLY(3  1=Y( U1+Y(21)+A.»Y(20) 

DO  56  K=12,  18,2 

PCLY13)  =  POLY(3)+'i.»YtK)+2.»Y(K+l) 

P0LY(3)=PULY(3)/30. 

WKITE  TAPE  8, I PROB , POLY, N, E I GENV , G.NREG, NBC, C, S IGMA, REPLTH, E IGENC, 

GC, ACCUR,HC,R,B 

GO  TO  (  15,26), IPROB 

WRITE  TAPE  0,X,Y 

CALL  INQUIR 

IF(SENSE  SWITCH  1)50, ^,9 

1F( POLYI 1) )50, 50,51 


C 

THIS  SECTION  COMPUTES  THE  ADPITII 

c 

PURPOSES 

50 

NSUB=0 

GC  TO  52 

51 

NSUB=3t(MREG-l)»A 

52 

WRITE  TAPE  8,  NSUB 

IF(NSUB)26,26,53 

53 

00  100  NFUNC=l,7 

NPERGP=26 

GC  TO  (  18,20,22,23,2'!,  1,2), NFUNC 

18 

DX=.04 

X(l  )=0.0 

GO  TO  39 

20 

GO  TU( 100,39),NREG 

22 

Xll)=-l.O 

DX=.08 

ARG=ROUTER 

IREG=NREG 

GO  TO  107 

23 

X(1)=0.0 

60 

NPERGP=NMAX/NP1 

IF(NPtRGP-26)59,59,58 

58 

NPERGP=26 

59 

DX=1.0/FLUAT(NPERGP-1) 

GC  TO  39 

2^ 

GO  TO  (  100,60)  ,.MREG 

1 

GO  ru  (100,3),NREG 

3 

1REG=1 

ARG  =  R 

X(l  )=-1.0 

DX=.D8 

GC  TO  10? 

2 

GO  TO  (  100, 'i)  ,NREG 

ti 

IREG=2 

10  7 

DO  108  K=1,NMAX 

119 


10!)    Y(K)=0.0 

i')    00    lOA    K  =  2,NPErtGP 
IC    X(Kl=X(K-l)+DX 

ir(NFUNC-2)'.2,A2,'.l 
'tl    DC     18    L=1,NP1 

TEMP=PHIL1L-1,A«G, IREG) 
DC  3a  K=1,NPERGP 
109  GC  TO  (31,32,33,31,32,36,36) ,NFUNC 

30  CONTINUE 
GO  TO  105 

'.2  DO  106  K=1,NPERGP 
GO  TO  109 
106  CONTINUE 
GO  TO  105 

31  ARG=X(K)»R 
IREG=l 
IF(NFUNC-1)57,57,34 

32  ARG=R+XIK)»REFLTH 
IREG=2 
lF(NFUNC-2)57,57,34 

57  YIKj^'i.^PI'PHILtO.ARG,  IREG) 
GO  TU  106 

33  Y(K)=Y(K)+P(L-l,X(K) )»TEMP 
GOTO  38 

S'i  NSUB  =  K<-(L-l)»NPERGP 

X(NSUB)=X(K) 

Y(NSUB)=PHIL(L-1 .ARG.IREG) 

GO  TO  38 
36  Y(K)=Y(K)+P(L-li X(K) )»TEMP 

GC  TO  38 
105  IFINFUNC-A)4't,'.3,'>3 

43  IF(NFUNC-6)A7,44,44 
47  NPTS=NPt»NPERGP 

NPLOTS=NPl 
GO  TO  45 

44  NPTS=NPCRGP 
NPLUTS^l 

45  YMIIM=0.0 
YMAX=0.0 

DO  46  K=1,NPTS 
YMIN=M1NI(YMIN,Y(K) ) 

46  YMAX=MAX1(YMAX,Y(K) ) 

WRITE  TAPE  8,X(1 ) , X ( NPERGP ) , YM IN , YMAX , I  SCALE , JSCALE , NPTS , NPLOTS , NO 
lOPY.NCARDS 
WRITE  TAPE  8,X,Y 
100  CONTINUE 
26  RETURN 
ENO 

BOP   PHIL  14 

FUNCTION  PHIL(L,X, IREG) 
THIS  SUBPROGRAM  COMPUTES  THE  LTH  MOMENT  OF  THE  ANGULAR  FLUX  AT  X. 


OIMCNSIGN  POLY  (6)  ,  E  IGEIMV  (  5  )  ,  G  (  1 1 ,  5  ) ,  GC  t  1 1 ,  5  )  ,D  (  2  I  t  SIGMA  (  2  )  ,E1GENC( 
15),BC(  U.bl  ,BI  10) 

CCMCON  PQLY,N,ElGENV,G,NREG.NBCiD,SIGMA,REFLTH,EIGENCfGC,ACCUR.8C, 
LK,B 

NCV2=(Ntl)/2 

PHIL=0.0 

GO  TO  (1,5),1REG 

1  DO  A  K=l,NGV2 
IF(K-1)2,2,3 

2  1  =  1 

GO  TO  4 

3  1  =  0 

1^    PHIL  =  PHIL  +  B(K)«GC(L+l,K)»FLaATI2»L+l)»C(L,SIGMA(l)»X/[  IGENCIK),  1  ) 
RETORN 

5  DO  6  1=1,2 

DO  6  K=1,N0V2 
NSU[S  =  K  +  N0V2*I 

6  PHIL=PH1L+B(NSUH)»G(L+1,K1»FL0AT(2«L+1)»Q(L,((-1.1»«(I-1) ) 'SIGMA {2 
l)«X/EIGENV(Kn«(-l.  )«»(L»1  I-l  )  ) 

R£TURl^l 
END 

POP   RGLTZMANN6  I'l 

;      THIS  PROGRAM  IS  THE  CONTROL  PROGRAM  FOR  PHASE  6 

DIMENSION  P0LY(6),E1GENV('>),G(  11,5),GC(11,5),C{2),SIGM«(2)  .EIGENCI 
l5  1,t)C(  11,5),B(  10) 

COMMON  POLY, N.EIGENV.G.NREG, NBC, C, SIGMA, REEL TH,e I GENC.GCACCUR.BC, 
IR.R 
REWIND  A 

READ  TAPE  I.NCASES 
DC  104  ICASE=1,NCASES 
TYPE  ^.ICASE 
\>    FORMAT!  16X1'.HPR0BLEM  NUMBER, T5) 
WRITE  OUTPUT  TAPE  6  ,  1 1 1 1 ,  ICA";E 

IIU  FORMAT  (2H 11,  15, 75H 

1 ) 

READ   TAPE  't ,  I  PRQB  ,  POL  Y  ,  N  ,E  I  f  ENV  ,  G.NREG  ,NBC,  C  ,  SIGMA,  RCFLTH.E  IGENC  , 
lGC,ACCUR,BC,R,a 
GO  TO  ( lOJ.lOO), IPRQB 

100  DO  101  K=l, U 

101  READ  INPUT  TAPE  5,6 
6  F0KMAT(80X) 

GO  TO  104 
101  CALL  OUTPUT 
104  CONTINUE 

REWIND  4 

CALL  EXIT 

STOP 

END 

HOP   OUTPUT  14 

SUHROUriNt  UUTPUI 
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C      THIS  SUBKQUriNE  PRliJTS  OUT  AN'D  PUNCHES  CARDS  FOR  ALL  PFRTINANT 
C        CALCULATED  DATA 

DIMENSION  P0LY(6),ElGENV(5),r(U,5),GCIll,5)  ,C12I  ,  S  I  GMM  2  )  ,  E  IGtNC  ( 

ir>),BC(  ii,5),n(  10) 

DIMENSION  X( 100) ,Y( 100) 
DIMENSION  YP(  in  ,  IP(  U) 
COMMON  P0LY,N,EIGENV,G,NREG,N'BC,C,SIGMA,REFLTH,E1GENC,GC,ACCUR,BC, 

IR.B 
COMMON  X,Y 
C      THIS  SECTION  PRINTS  AND  PUNC^^ES  THE  NORMAL  DATA 
NP1=N+1 
NCV2=(N<-n/2 
X(l)=SIGMA(l)«R 
WRITE  OUTPUT  TAPE  6,l,R.X(l) 

1  F0RMAT(3H0$0,20XilHTHE  COMPU-'ED  CRITICAL  RADIUS  I  S,F9.4, 1X3HCM./ 3H 
I  i     ,22X3II0R,F9.'.,  IX16HMEAN  F^EE  PATHS.) 

WRITE  OUTPUT  TAPE  6,121,P0LYn) 
121  FORMATOH  $  ,20X41HTHE  INTEGnATEO  INWARD  ANGULAR  FLUX  AT  THE/3H  $ 
1,22X19HVACUUM  INTERFACE  I S, 2 V IPE 10. 3 ) 
WRITE  OUTPUT  TAPE  6,2 

2  F0RMAT(3H-$-,  13X11HIN  THE  CWE:) 
WRITE  OUTPUT  TAPE  6,3 

3  FCRMAT(3H0$0,15Xl2HEIGeN  VALUES) 

WRITE  OUTPUT  TAPE  6,  A  ,  (  C  I  GENf"  (  K  )  ,  K=  I ,  NUV2  ) 

4  FDRMAT(3H  $  ,  17X  IPE  I  2  .5  ,  2XE  I  ■> .  5  ,  2XE  12  .  5  ,  2XE  1  2.  5  ) 
WRITE  OUTPUT  TAPE  6,5 

5  FCRMAT(3H0$0,  15Xr3HG(L,K)  MATRIX) 
DO  6  L=l,N0V2 

6  WRITE  OUTPUT  TAPE  6, 'i ,  (  GC  (  K  ,  L  )  ,K=  1  ,NP  1 ) 
GO  TO  (7, 10) ,NREG 

7  IF(NBC-A) 120, 10  ,120 
120  WRITE  OUTPUT  TAPE  6,8 

8  FORMAT!  3HOJ.Q,15X2!>HBOUNDARY  rONDlTION  MATRIX) 
DO  ')    L=l,N0\/2 

9  WRITE  OUTPUT  TAPE  6, 'i ,  (  BC  (  K  ,  I  )  ,  K=  1 ,  NPl  ) 

10  WRI  IE  OUTPUT  TAPE  6,11 

11  FORMAT (3H0$0, 15X23HN0RMALIZEC  COEFFICIENTS) 
WRITE  OUTPUT  TAPE  6, 4 , ( D ( K ) , K= 1 , N0V2 ) 

GO  TO  (  16, 12) ,NREG 

12  WRITE  OUTPUT  TAPE  6, 13 

13  F0RMAT(3H-$-,13X16HlN  THE  REFLECTOR) 
WRITE  OUTPUT  TAPE  6,3 

WRITE  OUTPUT  TAPE  6 , 4 , ( C I GENV ( K ) , K= 1 ,N0V2 ) 

WRITE  OUTPUT  TAPE  6,5 

DO  14  L=1,NQV2 

14  WRITE  OUTPUT  TAPE  6, 4 , ( G ( K , L ) , K= 1 ,NP 1 ) 
WRITE  OUTPUT  TAPE  6,8 

DO  15  L=1,N0V2 

15  WRITE  OUTPUT  TAPE  6 , 4 , ( BC ( K , L ) , K= 1, NP I ) 
N0V  =  N0\/2-i-l 

NSUB=3«N0V2 

WRITE  OUTPUT  TAPE  6, 11 
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WRITE  OUTPUT  TAPE  6, 4 , ( B ( K ) , K=NQV ,NSUB ) 

16  WRITE  OUTPUT  TAPE  6, 17 

17  FCRI'!AT(3H-$-,13X4BH0ATA  FOR  ANGULAR  FLUX  PLOT  AT  THE  OUTER  BOUNDAR 
lYI 

WRITE  OUTPUT  TAPE  6,18 

18  FCRKATl  3HQ$0,17X,3(2HMU,'.X7HPHI  IMU)  ,6X)  ) 
READ  TAPE  'i.X.Y 

DO  20  K=l,7 

20  WRITE  OUTPUT  TAPE  6, 21 , ( X ( L) , Y ( L ) ,L=K , 21 ,7 ) 

21  F0RMAT(3H  $  ,15X,3(0PF't.l,2X,lPEl0.3,3X)  1 
READ  TAPE  'i ,     NTIMES 

IF(NTIMES)  106, 106, 105 
105  CALL  INQUIR 

IFISENSE  SWITCH  1)108,109 

108  P0LY(2)=0.0 

THIS  SECTION  PRINTS  THE  EXTRA  OATA  USED  FOR  TESTING  PURPOSES 

109  DO  50  M=l,7 
IF(NTIMES-3)61,61,73 

61  IF(M-2)7l,72,71 

71  lF(M-5)73,72,72 

72  IF(P0LY(2) )50,50, 110 

110  00  63  K=l,2 

63  READ  INPUT  TAPE  5,64 

64  rORMAT(aOX) 
GO  TO  50 

73  READ  TAPE  4 , XM IN , XMAX , YM I N, YKAX ,  I  SCALE , JSCALE , NPTS.NPLOTS, NCOPY , NC 
lARDS 

READ  TAPE  4,X,Y 

GO  TO  (80,24, 38, 81,40,96,99) ,M 
80  WRITE  OUTPUT  TAPE  6,90 

90  FORMAT! 1H1,30X73HTOTAL  FLUX  DISTRIBUTION  IN  THE  CORE  AS  A  FUNCTION 
1  OF  THE  RADIAL  DIMENSION///) 

WRITE  OUTPUT  TAPE  6,22 

22  FORMAT! 15X1HR,16X6HPHI(R)//) 
00  23  K=1,NPTS 

23  WRITE  OUTPUT  TAPE  6, 28 , X ! K ) , v ( K ) 
GO  TO  74 

24  WRITE  OUTPUT  TAPE  6,91 

91  FORMAT! 1H1,26X78HT0TAL  FLUX  PISTRIBUTION  IN  THE  REFLECTOR  AS  A  FUN 
ICTION  OF  THE  RADIAL  DIMENSION///) 

WRITE  OUTPUT  TAPE  6,22 

00  36  K=1,NPTS 
36  WRITE  OUTPUT  TAPE  6 , 28 , X ( K ) , Y ! K ) 

GO  TO  74 
38  WRITE  OUTPUT  TAPE  6,92 

92  FORMAT! IH1,37X56HANGULAR  FLUX  AT  THE  VACUUM  INTERFACE  AS  A  FUNCTIO 
IN  OF  MU///) 

WRITE  OUTPUT  TAPE  6,25 
2  5  FORMAT ! 13X2HMU, I 2X14HPHI ! ROUTER , MU) // ) 

DO  26  K=l,NPrS 
26  WRITE  OUTPUT  TAPE  6 , 28 , X ! K ) , Y ( K ) 
26  FORMAT! 11X0PF6.3, 12X1PEI0.3) 
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GC  TO  T< 
81  DO  J-V  K=1,NP1 
3^1  IP(K)  =  K-1 

WRITE  OUTPUT  TAPE  6,93 

93  FORMAT! 1H1,30X69H1ND1VI0UAL  SPATIAL  MOMENT  OISTKIBUTICN  IN  THE  COR 
IE  AS  A  FUNCTION  OF  K///) 

95  WRITE  OUTPUT  TAPE  6, 33 , ( I P ( K ) , K= 1 ,NP 1 ) 
33  FDRKAT(7X6HR  ,L  = , IX , 1 1 ( I  3, 8X ) // ) 

NPERGP^NPTS/NPLOTS 
DO  31  K=1,NPERGP 
DO  30  L=1,NP1 
NSU8=K+(L-11«NPERGP 

30  YP( L)=YINSUB) 

31  WRITE  OUTPUT  TAPE  6, 32 , X ( K ) , ( YP ( L I , L=l ,NP1 ) 

32  F0RMATt5X,F5.2,ll( IXIPEIO.3) ) 
GO  TO  TA 

^.0  WRITE  OUTPUT  TAPE  6,94 

94  F0RMAT(1H1,29X74HIN0IV1DUAL  SPATIAL  MOMENT  DISTRIBUTION  IN  THE  REF 
HECTOR  AS  A  FUNCTION  OF  R/// 1 

GO  TU  95 

96  WRITE  OUTPUT  TAPE  6,97 

97  FORMAT! 1H1,30X72HANGULAK  FLUX  DISTRIBUTION  AT  THE  INTERFACI&L  BOUN 
IDARY  FROM  THE  CORE  SIDE///) 

101  WRITE  UUTPUT  TAPE  6,41 
41  FORMAT! 13X2HMU,15X9HPHIIR,MU)//) 
DO  111  K=1,NPTS 
111  WRITE  OUTPUT  TAPE  6 , 28 , X ! K ) , Y I K I 
GD  TO  74 
99  WRITE  OUTPUT  TAPE  6,100 
100  FORMAT! 1H1,28X77HANGULAR  FLUX  DISTRIBUTION  AT  THE  INTERFACIAL  BOUN 
IDARY  FROM  THE  REFLECTOR  SIDE///) 
GC  TU  101 
74  IF1P0LY(2)  )ia7,50, 107 

THIS  SECTION  PRINTS  THE  GRAPHS  IF  THEY  ARE  CALLED  FOR 
10  7  CALL  PL0T!X,Y,XMIN,XMAX,YMIN,YMAX,ISCALE,JSCALE,NPTS,NPLOTS,NCQPY, 

INCAROS) 
50  CONTINUE 
106  RETURN 
END 

GOP   PLOT  14 

THIS  SUBPROGRAM  PLOTS  THE  GRAPHS 

SURaUUTINLPLUTIX,Y,XMIN,XMAX,YMIN,YMAX,LX,LY,NPT,NPLUT,NCOPY,NCD) 
D1MENSI0NX!1),Y( 1),SX( 10),TITLEI8),TABI4),LI 135) ,NCH(in),MOPI30) 
NC0=NCD+1 
GOTO! 1,00,82,80) ,NCD 

80  READ  INPUT rAPE5,81, (TI TLE !  I T ) , I T= 1 , 8 ) 

81  F0RMATI8A10) 
IFI4-NCD)1,82, 1 

8  2  READ) NPUTI APES, 83, !MOPI  1) , 1=1,30) , I NCH ( I ) , I  =  1 ,  10) , ( TAR  I  I ) ,1  =  1,4) , 

IND,NP,NM,NB 
83  FORMAT!40A1,4A9,4A1) 
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1 

0=    2.3025B5 

IF(NPL0T-1 )2,2,3 

2 

NPN=NPT 

GO    TO    4 

i 

NPN=NPT/NPLOT 

4 

IF(LX)5,5,7 

5 

CX=120./(XMAX-XMIN) 

SX(  1)  =  XMIN 

SXI6)=XMAX 

Z=XMIN 

D06K=2,5 

Z=(XMAX-XMIN)/5.+Z 

6 

SX(K)=Z 

00    TU    11 

7 

CX=120/LX 

NX=1L0G(XMIN) )/0 

LLX  =  LXH 

D010K=1,LLX 

10 

SXtK)=10.»»(NXtK-l) 

11 

IF(LY)12,12,13 

12 

CY=50./(YMAX-YMIN) 

GO    TO    16 

13 

CY=50/LY 

KY=CY 

1 

NY=(LQG(YMIN) ) /Q 

'                      16 

IF(LX)  17,17,21 

17 

00201=1, NPT 

IF(XMIN) 18,19,19 

18 

M=CX«X( I )+.5-CX»XMIN 

GCTCJ20 

19 

M-CX«X(  I  1  +  .5 

20 

X(I )=M 

GCT023 

21 

00221=1, NPT 

M=CX»(LOG(X( n/XMIN)/Q)+.5 

22 

X(I  )  =  M 

2  1 

IF(LY)2'.,2A,28 

2A 

00271=1, NPT 

IF(YMIN)25,26,26 

25 

M=CV«Y( I  )t.5-CY»YMIN 

GO    TO    27 

26 

M  =  CY»Yl  I  1  +  .5 

27 

Yd  )  =  M 

GO     ro    30 

28 

DO    291=1, NPT 

M=CY»(LOG(Y( I)/YMIN)/0)*.5 

29 

Y( I )=M 

30 

DO    79    NN=1,NC0PY 

M=l 

T1=^0. 

LYY=LY 

TT=50. 

aw 
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WKITEOUT PUTT  APES, 31, ( T ITLC ( I T ) , I T  =  l , 8 ) 

31  FCRHAT( 1H126X8A10) 
0C61  KK=1,51 

N=l 

NNN=NPN 
JED=l 
T=51-KK 
D032J=1, 132 

32  L(J)=NB 
L(132)=N0 
L( 12)=N0 
IF(LY)33, 33,38 

33  L( 12)=NP 
IF(T-TT)3't,3'i,'.0 

3'.  SCALE=T/CY 
L(132)=NP 
IF(YMIN)35,36,36 

35  SCALE=T/CY+YMIN 

36  N  =  0 
TT=TT-5. 

IF(  T)37,37,'.0 

37  SCALE=YMIN 
GCT040 

38  SS=KY»LYY 
IF(T-SS)39,39,A0 

39  SCALE=10.»«(NY+LYY) 
N=0 

LYY=LYY-1 

L{12)=NP 

L(132)=NP 
'tO  IF(50.-T)A3,4'.,',3 
ifi    IF{T)50,'.A,50 
',',    D045J=13,  132 
^5  L( J)=NM 

IF(LX)46,^6,'.8 

46  D047J=12, 132, 12 

47  L(J )=NP 
IFCiO.-TjbC'il.fiO 

41  L(132)=ND 
G0T05O 

48  KX=120/LX 

DO  49  J=12,132,KX 

49  L(J)=NP 
IF(50.-T)'jO,42,50 

42  L(132)=NU 

50  Da53LM=l,NPL0T 
D(j521  =  JCO,N.NN 
IFlYd  )-T)52,51,52 

51  J  =  X( I  1  +  12. 

1F( J-12)52,85,84 

84  IF( J-132)85,85,52 

85  L(J)=NCHtLMI 
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52    CONTINUE 

JED=NNN+1 

NNN=NNMtNPN 
5  3  CCNTINue 

IF( 11-1)56,5^,56 

54  IF(  10.-1)55,56,56 

55  L(2)=M0PIK1) 
Ml=M+l 
Tl=Tl-l. 

56  IF(N)57,57,59 

5  7  WRI TE0UTPUTTAPE6, 58,L I  2), SCALE, (LlJltJ'l 1.1321 
58  FORMATdH  A  I ,  E8.  2  ,  122A  I ) 

G0TU61 
54  WRITE0UTPUTTAPE6,60, (L(J),J=1,132I 

60  FORMATt 132A1) 

61  CCNTINUE 
IF(LX)6A,62,64 

62  WRITE0UTPUTTAPE6,63,(SX(K),K=1,6) 

63  FCRMAT(7XE9.3,l5XE9.3,l5XCy."?,l5XE9.3,l5XE9.3,UXE9.3) 
GG  ru  77 

64  GO  TO  (65, 67, 69, 71, 73, 75). LX 

65  WRITE  aUTPUTTAPE6,66, (SX(K),K=1,LLX) 

66  F0RMAT(7XE9.3,107XE9.3) 
GO  TO  77 

67  WRITE  0UTPUTTAPE6,68,(SX(K),K=1,LLX) 

68  FGRMAT(7XE9.3,52XE9.3,'.6XE9.?) 
GO  TO  77 

69  WRITE  0UTPUTTAPE6,70, ISXIK),K=1,LLX) 

70  FQRMAT(7Xt9.3,3lXE9.3,31XE9.T,2  7X69.3) 
GO  TO  77 

71  WRITE  0UTPUTTAPE6,72,(SX(K),K=l,LLX) 

72  FORMAT! 7XE9. 3, 21 XE9. 3, 2 1XE9. 1, 2 1 XE9. 3, I 7XE9. 31 
GO  TO  77 

73  WRITE  0UTPUTTAPE6,74,(SX(K),«=1,LLX) 

74  F0RMAT(7XC9.3,15XE9.3, 15XE9.', 15XE9. 3. 15XC9. 3, 1 1 XE9. 3 ) 
CO  TO  77 

75  WRITE  0UTPUTTAPE6,76,  (SX(K),''=1,LLX) 

76  F0RMAT(7XE9.3,UXE9.3,  1 1XE9.  '  ,  1  1  XC9.  3,  1 IXE9.  3, 1  lXt9.  3  ,  7XE9.3  I 
7  7  WRIIEaUTPUTTAPE6,78, ( TAB ( 1 ) , '= I , 4 ) 

78  F0RMAT(48X4A9) 

79  CONTINUE 
RETURN 
END 
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A  study  was  made  of  tho   convergence  of  the   spl\erical   harmonics  approxi- 
mations of  both  even  and  odd  orders  to  the  monoeneruetic  Boltzmann  neutron 
transport  equation   for  a  nuclear  reactor.     The  theory  and  computer  programs 
required  for  determining  the   critical   radii   of  bare   and   reflected  spherical 
reactors  were  developed.     Considerable  care  was  taken  in  the  development  of 
the  tlieory  for  the  even  order  approximations   in  order  to  obtain   consistent 
results.      Pour  types  of  boundary  conditions  were  employed  at  the  vacuum  inter- 
face of  the   reactor  system,   namely;      Marshak's,   Mark's,   replacement   of  the 
vacuum  with  an   infinite   black  medium,    and  a   set   of  boundary  conditions 
derived   from  a  variational    formulation   of  tlie   Boltzmann  equation. 

The  nature  of  the  convergence  of  the  spherical  iiarmonics  approximations   for 
a  bare   spherical   reactor  was   found  to  be  dependent   on  the  boundary  conditions 
applied.     The   critical   radii  wliich  resulted  from  using  Marshak's  boundary  con- 
ditions converged  asymptotically  from  above   along  separate  paths   for  the  even 
and  odd  order  approximations.      For  Mark's   and  the   infinite  black  reflector 
boundary  conditions  tlie  even  order  approximations   converged  asymptotically  from 
below  while  the  odd  order  approximations   converged  asymptotically   from  above, 
i.e., they  counterconverged  to  the  exact   result.     No  conclusions  could  be  drawn 
concerning  the  nature   of  tlie  convergence  wlien  variational  boundary  conditions 
were  employed.     Whenever  one  or  more  of  three  special   conditions  were  satisfied 
in  a  given   reflected  reactor  system,   the  computed  critical   radii  were   found  to 
counterconverge  to  the  exact   result.     The  odd  order  approximations  were   found 
to  be  superior  to  the  even  order  approximations   for  all  but  the   lowest  orders 
of  approximations.     The   1\   approximation   almost   always  yielded  results  better 
than  tliose   given  by  a  P     approximation. 


When  the   critical   radii   resulting  from  the   splierical   harmonics  approxima- 
tions were   found  to   counterconverge  to  tlie  exact   result   a  particular  weighted 
average  of  the  even  and  odd  order  results  was   conputed.     This  weighted  average 
was  found  to  yield  a  considerably  more  accurate  estimate  of  the  critical   radius 
than  either  of  the   individual  results  when  approximations  of  order  three  or 
greater  wore  utilized.     The  odd  order  approximations  using  variational  bound- 
ary conditions  were   found  to  yield  the  best   results   for  a  given  order  of 
approximation.      Althougii  many  more  cases  need  to  be  examined,   the   results 
obtained  here  strongly  indicate  that  tlie  even  order  approximations  are  of  con- 
siderable value. 


